1. Introduction

Photonic crystal fibers (PCFs) [1, 2] are novel highly efficient components applied in nonlinear optics [3, 4], ultrafast science [5], optical metrology [6], nonlinear spectroscopy [7] and microscopy [8], biomedical optics [9], and optical sensing [10]. Unlike standard fibers, waveguiding in PCFs is often a result of interference of multiple light waves scattered and reflected by refractive index discontinuities in an air-silica microstructured cladding. This nature of the guided modes in PCFs makes them interesting objects of optical physics and offers convenient and practical ways to tailor the field intensity profiles and the dispersion properties of the light waves and the laser pulses propagating through PCFs [11].

Although some important properties of the guided modes in PCFs can be qualitatively understood with the use of analytical or semi-analytical approaches [12], their quantitative analysis generally requires rigorous numerical methods. Several powerful fully vectorial techniques have been recently adapted to simulations of light propagation in various classes of PCFs. The most widely used and successful approaches include the plane wave expansion [13, 14], the polar Fourier decomposition method (FDM) [15], the finite element expansion (FEM) [16, 17], the expansion in localized functions [18], the multipole technique [19, 20], the source model technique [21], and the FDTD method [22]. A combination of these approaches has provided important insights into fundamental aspects of the guided mode field intensity profiles and their relation to the fiber symmetry, the dispersion parameters, and the various types of losses.

 

Fig. 1. SEM image of the nonlinear PCF.

Download Full Size | PDF

In this work, we present a comparison of four of the above mentioned fully vectorial methods by calculating the spectral dependence of the phase and the group effective indices, the birefringence, and the group velocity dispersion for three spatial modes of the lowest order guided in a PCF with large holes of noncircular shape. The calculations are performed for the real fiber structure shown in Fig. 1. It is very important to reproduce the fiber geometry in an unambiguous and accurate way in order to avoid possible discrepancies between different methods related only to differences in the assumed input fiber geometry. To this aim, we process the SEM image of the fiber cross section by modifying its histogram and applying thresholding and binarization procedures to detect the edges of the cladding holes. As a result, we obtain a binary image of the fiber cross section reflecting the real shape and the location of each hole with the accuracy of about 40 nm. In the last step, the edge of each hole is smoothed by applying a cubic spline approximation. This last fiber geometry defined by splines is used as input data in all the calculation methods. Still, calculation errors related to the precision of the specific expansion technique used to represent the input fiber geometry are possible. The comparison of the simulation results obtained with the four different rigorous methods provides therefore useful information about the practically achievable accuracy in modeling the propagation characteristics of real PCF structures.

2. Rigorous methods

Four fully vectorial methods, i.e., the finite element method (FEM), the localized function method (LFM), the plane wave method (PWM), and the source model technique (SMT) have been used to model the characteristics of guided modes in the PCF shown in Fig. 1. Earlier experiments [23] have demonstrated that this type of PCF can serve as an efficient novel light source utilizing nonlinear optical spectral transformation of ultrashort laser pulses. The dimensions 4.20×3.39 µm of the elliptically deformed (but not precisely elliptical) core of this fiber were determined using the procedure described in section 3. One and the same PCF geometry represented by a spline approximation is used as input data in all numerical methods. Furthermore, in all methods, we take into account the chromatic dispersion of the silica glass represented by Sellmeier equation [24].

Table 1. Results of the convergence test carried out for the ${\mathit{\text{LP}}}_{\mathit{01}}^{\mathrm{x}}$ mode at λ=1.55 µm.

View Table | View all tables in this article

2.1. The finite element method

The first numerical approach employs a finite element representation of the analyzed structure (FEM) [23]. In this technique, a full vector mode solver based on a hybrid edge/nodal finite element method (FEM) [17, 25] is used to determine the propagation constants and the electric field distributions of the guided modes by solving the following generalized eigenvalue problem (GEP):

Here, the eigenvector is the electric field vector in the PCF cross-section, $\overrightarrow{E}$ =[$\overrightarrow{E}$ _⊥,E _z]. The eigenvalue is the mode propagation constant, β, and k ₀=ω/c is the wave number. Throughout this work we assume that the z and time dependence of all fields is exp[j(ωt-βz)]. In order to find the eigenpairs of the above equation, we use an Arnoldi method as an eigensolver and a unsymmetric multifrontal method as a solver of a linear system of equations. Using the spline representation of the fiber cross section as input data, we generate the mesh for finite element calculations (Fig. 2). The mesh is composed of 630 K triangular elements and divided into two regions with different average size of the elementary triangle equal respectively to 19 nm and 100 nm. By performing a convergence test for the effective index (n=βc/ω) of the fundamental mode at λ=1.55 µm, we find that for the mesh composed of 630 000 the calculated value of n is stable on the sixth decimal place, see Table 1.

 

Fig. 2. Mesh used in finite element method (a) and super cell used in plane-wave method (b) with indicated co-ordinate system.

Download Full Size | PDF

2.2. The plane wave method

The second rigorous approach used in the simulation exercise is based on a plane wave method (PWM) [13, 14]. The details of our implementation are summarized in [31]. Briefly, we are solving the eigenequation for the transverse magnetic field, which unlike the electric field is continuous at the refractive index discontinuities and may be represented more accurately using the same number of continuous basis functions:

Equation (2) retains the form of an eigenequation also for dispersive materials, because we place the propagation constant in the position of the eigenvalue. This comes in contrast with the other popular PWM formulation [14] where the frequency is the eigenvalue, the treatment of dispersion is not direct however, the eigenequation is Hermitian. We use Arpack to solve the resulting algebraic eigenproblem. This requires providing a sequence of matrix-vector products that we evaluate in a computationally efficient way using FFT.

As the microstructure of the fiber contains no periodic elements, there is no obvious preference for choosing a specific supercell shape or size. On the other hand, the trade-off between its size and the accuracy of representation with a given number of plane-waves, indicates that using a large supercell size may be counterproductive. We decided to apply a rectangular supercell of dimensions 15.00×13.86 µm, comprising the central part of the fiber cross-section [Fig. 2(b)]. The supercell of such size contains only very small areas of the solid part of the cladding located outside the second layer of holes. This is not a problem because due to the large diameters of the cladding holes, the amplitudes of the guided modes drop abruptly outside the core (see Fig. 5). Therefore, the chosen supercell is sufficiently large for accurate simulations of the guided modes. However, the periodic boundary conditions imply that the PWM is accurate to the limit of the radiation losses, which are neglected in this case. The discontinuities of the index profile resulting from connecting the supercells into an infinite periodical structure does not compromise the accuracy of plane-wave representation, since the same type of discontinuities exist in the center of the PCF, where the modal amplitudes are much higher.

The Fourier decomposition of the permittivity profile is calculated numerically based on the intermediate high resolution raster representation consisting of square pixels with the size of Δl=20 nm. Apparently, the intermediate discrete representation of the PCF cross-section is the major limiting factor for the accuracy of our PWM calculations. We use 2¹⁴ plane waves for each of the two Cartesian field components.

The convergence analysis given in Table 1 shows that for such number of basis functions the effective index n=β/k ₀ is stable to the sixth decimal place at λ=1.55 µm. On the other hand, the order of error introduced by the intermediate raster representation may be roughly estimated as Δn=Δl·dn/dλ. For λ=1.55 µm, this error is equal to 10^-3 indicating that the PCF under study is highly sensitive to distortions or manufacturing precision.

Furthermore, convergence problems are also clearly visible in the results from Table 1, like it is often the case with the PWM [14, 26, 27]. Here, the slow convergence manifests itself through oscillations that may be attributed to the presence of high frequency components in the spatial spectrum of the structure. These components result from the presence of fine details in the considered PCF. They are included in the calculations only after a sufficiently large number of plane-wave basis functions are taken into account.

2.3. The localized function method

The third numerical approach is based on the modified fully vectorial localized function method (LFM), following the generic algorithm developed by Monro et al. [18]. It solves the vectorial wave problem for the electric field, E⃗=[E⃗_⊥,E _z], in the PCF cross-section. The spline representation of the realistic n ²(x, y) distribution of the PCF structure is accurately approximated with 150×150 trigonometric expansion terms and 25×25 Hermite-Gaussian polynomials. The E _x and E _y components of the electric field are represented as series expansions in 25×25 Hermite-Gaussian polynomials. The substitution of the series expansions for E _x, E _y, and n ²(x, y) into the wave equations for the transverse field components reduces the vectorial wave equation to an eigenvalue problem and allows to determine the propagation constants and the transverse field profiles of the guided modes. The convergence test shows that the stability of the effective index in the fourth decimal place is achieved at λ=1.55 µm, (Table 1).

2.4. The source model technique

A source model technique (SMT) has been recently developed for the analysis of PCFs [21, 28]. Like the FEM, the LFM, and the PWM, the SMT is a fully vectorial method that can be used to determine the modes of a cylindrical structure of piecewise-homogeneous cross section.

In the SMT, the fields in each homogeneous region of the cross section are approximated by the fields due to a linear combination of elementary sources placed outside of it, on curves conformal with the region boundary. The sources radiate in a homogeneous medium with the same material parameters as those of the region they enclose in the PCF cross section. Their fields, therefore, have well known analytic expressions. The amplitudes of the elementary sources are determined so as to satisfy the continuity conditions across the media interfaces at a set of testing points. Usually, more testing points than sources are used, and the solution is sought in the least squares sense. The elementary sources used are electric and magnetic current filaments, carrying longitudinally varying currents, that vary with the z-coordinate as exp(-jβz). For each pair of (ω,β) under consideration, the vector that yields the smallest least squares error (normalized to the root-mean-square value of the modal fields) is found via the solution of a GEP, in a way similar to that described in [29]. Once a low error solution is found for some (ω,β) pair, the current amplitudes are inserted into the linear combinations that approximate the fields, yielding the mode field patterns and any other parameter of interest.

It should be noted that the GEP is quite different from the GEPs that occur in the FEM, the LFM, and the PWM. In the latter methods, the eigenvalue is β (or ω in some formulations), whereas in the SMT, the eigenvalue is the error in the continuity conditions. The matrices in the FEM and PWM do not depend on ω or β (whichever is the eigenvalue). In the SMT, because the basis functions used inherently satisfy Maxwell’s equations throughout the PCF cross-section, the matrices naturally depend on both ω and β. The GEP must be solved a number of times, at different (ω,β) pairs chosen by an optimization algorithm, until a pair which yields a small enough error is found. Although this process appears to be a drawback of the SMT, it has some advantages. First, a reliable indication of the quality of the solution is readily available, namely the eigenvalue. Second, the representation of the field in terms of the SMT basis functions tends to be more compact than in the other methods, so each solution of the GEP requires less computational resources.

The SMT has much in common with the Multipole Method [19, 20], in which the basis functions employed, multipole fields of increasing orders, also satisfy Maxwell’s equations. Although this representation is considered very compact, it is directly applicable only to cross sections composed of circular homogeneous inclusions in a homogeneous medium. In the SMT, the homogeneous regions in the cross section can be arbitrarily shaped. In the context of PCFs, however, only circular [21, 28] and elliptical [28] air holes have been studied. The application of the method to the analysis of a realistic fiber, with its irregularly shaped air holes, requires the development of a systematic way to locate the sources around each homogeneous region.

It is important to be able to easily modify the number of sources and testing points, particularly for convergence tests. This can be done most conveniently by use of a parametric representation of the curves on which the sources and the testing points are to be distributed. Since the boundaries are approximated by splines, we already have a parametric representation of the air-hole boundaries, on which the testing points are distributed. A similar spline representation of the curve on which the sources are distributed is derived there from, by fitting a spline through the tips of normal vectors pointing outwards from each homogeneous region [see Fig. 3(a)]. The testing points and sources are distributed homogeneously on their respective curves, meaning that adjacent testing points and sources are separated by a constant arc length.

The fields of all the modes under consideration are strongly confined, so we only model the first ring of air holes, thus reducing the computational burden. At λ=1.55 µm, the addition of an outer ring alters the phase refractive index by approximately 3×10^-7, which is considered small enough to justify the approximation. When calculating confinement losses this approximation cannot be made, naturally, so we model the entire structure. In all the calculations except the convergence test, we use 50 electric current filaments and 50 magnetic current filaments placed around each air hole to approximate the fields there. Similarly, we place 50 sources of each type in each region originally occupied by an air hole to approximate the fields in the dielectric. The continuity conditions are enforced at 65 testing points per hole. The distance of the sources from the boundary is defined relative to the minimum radius of curvature of each hole, denoted by $r_{\text{min}}^{\left(i\right)}$ for the i-th hole. The sources used to approximate the fields in the air holes are placed at a distance of 1.2×$r_{\text{min}}^{\left(i\right)}$ i from the boundary, and those used to approximate the fields in the dielectric are placed at a distance of 0.65×$r_{\text{min}}^{\left(i\right)}$ from the boundary. The configuration of sources and testing points is shown in Fig. 3(b). As shown by a convergence test carried out at λ=1.55 µm, the applied number of sources and testing points assures a solution for effective index that is stable up to the sixth decimal place.

 

Fig. 3. Sources placed in a homogeneous air region that are used to approximate the fields in an air-hole and normal vectors used to construct the curve on which the sources are placed (a). Location of sources and testing points are shown in (b). Testing points are marked by x’s, sources used to approximate the fields in the air-holes are marked by hollow circles, while sources used to approximate the fields in the dielectric are marked by solid circles.

Download Full Size | PDF

3. An approximate equivalent fiber method

It has recently been shown in [32] that the propagation characteristics of a nonlinear PCF with residual core ellipticity can be effectively modeled using a simplified method, in which the real structure of the PCF is replaced by an equivalent step index waveguide composed of an elliptical glass core surrounded by an air cladding. Such an approach is used for the first time in [30] to model the spectral dependence of the phase and the group modal birefringence of the fundamental mode of a nonlinear PCF with a residual core ellipticity. The equivalent fiber method (EFM) seems to be intuitively well justified as long as fused silica bridges supporting the PCF core can be considered as small deviations from its perfect elliptical shape.

In this work, we verify the applicability of this approach to modeling the phase and the group effective indices, the birefringence, and the dispersion properties of the lowest order spatial modes supported by the photonic crystal fiber with the elliptically deformed core shown in Fig. 1. The first step in this method requires determining the dimensions of the equivalent fiber core. For this purpose, we use as objective criteria the overlap coefficient between the core of the equivalent elliptical fiber and the real fiber. This coefficient is calculated as a function of the sizes of the minor and the major axes of the equivalent elliptical core. As shown in Fig. 4, the highest overlap coefficient is obtained for a=3.39 µm and b=4.2 µm, which are the core dimensions used in the numerical analysis.

There exist several methods, including analytical [34], semi-analytical [33], and fully numerical [25] approaches that allow calculating the propagation constants of the conventional elliptical core fiber. In order to minimize the calculation errors and to compare the discrepancies related only to structural differences between the equivalent and the real fiber, we apply the FEM with a high mesh density (630 000 elements) to model the propagation characteristics of three spatial modes of the lowest order supported by the equivalent elliptical core fiber.

 

Fig. 4. Overlap coefficient vs. minor and major axes of the equivalent elliptical core (a) and the optimal elliptical core superimposed on the cross-section the considered PCF (b).

Download Full Size | PDF

As long as the structural differences between the real and the equivalent fiber can be considered as small deviations from the perfectly elliptical shape of the core leading to small perturbations of fiber modes, the correspondence between the modes of the equivalent fiber and the investigated PCF is more than just a formal convention. As a matter of fact it rests upon physical arguments that can be conveniently and instructively expressed in terms of the reciprocity theorem adapted to optical waveguides [34]. Our numerical results on the amplitude distributions, the birefringence, and the dispersion properties of the considered PCF, presented in the next section, provide helpful illustrations to these arguments. Therefore, we proceed with a more detailed discussion of the mathematical aspects behind the assignment of the PCF modes.

In its scalar version, the reciprocity relation for the field profiles Ψ and $\tilde{\Psi }$ , and the propagation constants β and β̃ of two waveguides with refractive index profiles n=n(x,y) and ñ=ñ(x,y), respectively, is written as [34]:

where the integration is performed over the entire plane of the fiber cross section.

Suppose that an equivalent elliptical core fiber with refractive index profile ñ(x,y) supports a family of guided modes $\tilde{\Psi }$ _j and radiation modes $\tilde{\Psi }$ _j (q) (q is the parameter of radiation modes). These modes are defined as solutions of the relevant wave equation and include the modes belonging to the LP family. We can then use, following [34], the known modes Ψ _j as a basis to expand the unknown modes Ψ _k in the real PCF with a refractive index profile n(x,y):

where a _j and a _j (q) are the unknown expansion coefficients.

As can be seen from Eq. (3), a small deviation of the refractive index profile n ²-ñ², gives rise to small perturbations of the guided wave solutions, allowing us to apply an iterative procedure to find the expansion coefficients in Eq. (4). In the zero order approximation, we set

Substituting Eq. (4) as zero order solutions into Eq. (4) and using Eq. (3), we arrive at

The above equations provide a convenient framework for the qualitative analysis of the results presented in the next sections.

4. Calculation results and discussion

The transverse distributions of the field intensity with electric field polarization maps calculated by the FEM for the low-order modes of the PCF are presented in Figs. 5(a)–5(c). The mode profiles obtained by the other methods are found to be in excellent agreement with those displayed in Fig. 5. The mode assignment in Fig. 5 and throughout the paper follows the classification adopted for conventional elliptical core fibers with low core ellipticity (0.5<a/b<1) [24]. We classify the fundamental mode with a bell-shaped field intensity profile reaching its maximum at the center of the PCF core [Fig. 5(a)] as the doublet of x- and y-linearly polarized LP ₀₁ modes, corresponding respectively to the ${\mathit{\text{HE}}}_{\mathit{11}}^e$ and ${\mathit{\text{HE}}}_{\mathit{11}}^o$ hybrid modes. In Fig. 5, we also present the field distributions for the modes with two-lobe spatial field intensity profiles, which are a pair of even LP ₁₁ modes (${\mathit{\text{LP}}}_{\mathit{11}}^{\mathit{\text{ey}}}$ and ${\mathit{\text{LP}}}_{\mathit{11}}^{\mathit{\text{ex}}}$ ) arising due to transformation of ${\mathit{\text{HE}}}_{\mathit{21}}^e$ and ${\mathit{\text{HE}}}_{\mathit{21}}^o$ hybrid modes, and a pair of odd LP ₁₁ modes (${\mathit{\text{LP}}}_{\mathit{11}}^{\mathit{\text{ox}}}$ and ${\mathit{\text{LP}}}_{\mathit{11}}^{\mathit{\text{oy}}}$ ) arising due to transition of TE ₀₁ and TM ₀₁ hybrid modes. The PCF modes differ slightly from their counterparts in a conventional elliptical core fiber.

 

Fig. 5. Field intensity profiles with polarization maps calculated with FEM for x-polarized LP ₀₁ (a), ${\mathit{\text{LP}}}_{\mathit{11}}^e$ (b), and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ (c) spatial modes.

Download Full Size | PDF

In particular, the removal of two-plane symmetry makes the amplitude distribution in the PCF modes more complicated and diversified compared to the modes of equivalent elliptical core fiber.

The numerical results presented later in this section demonstrate that the confinement losses for the considered PCF are low. This allows us to neglect in the expansion of Eq. (4) the coupling of the unperturbed guided modes $\tilde{\Psi }$ _j to the radiation modes $\tilde{\Psi }$ _j (q). The irregularities of the field intensity profiles of the PCF modes presented in Fig. 5 can therefore be viewed as a result of the coupling of the unperturbed LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes of the equivalent fiber to the higher order modes of the same fiber. These perturbations are dominated by the polar-angle-dependent high order modes that maximize, as visualized by the integrals in Eq. (6), the overlap with the exponentially decaying profiles of the unperturbed LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes and the high-index bridges supporting the fiber core. The spatial asymmetry of these bridges, represented by the difference n ²-ñ² in Eq. (6), is thus mapped on the asymmetry of the PCF modes shown in Figs. 5(a–c), leading to interesting birefringence properties of the PCF fiber that are analyzed in detail in the following paragraphs.

We perform simulations of several parameters of the considered PCF such as the propagation constants β, the phase and group indices defined respectively as n=βc/ω and n _g =c∂β/∂ω, the phase and the group modal birefringence, respectively Δn and Δn _g , the group velocity dispersion (GVD) parameter k ₂=d ² β/dω ², and the confinement losses. The simulations are carried out for the three lowest order spatial modes (LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ ) using four rigorous approaches outlined in the previous section, i.e., the FEM, the PWM, the LFM, the SMT, and the approximate equivalent fiber method (EFM).

The differences between n and n _g calculated using the different methods are lower than 4×10^-3, therefore, we show in Fig. 6 only the FEM results in order to avoid obscuring the plots. In Table 2, we gather the numerical results for all the simulated parameters for λ=0.4 µm and 1.55 µm.

The spline representation of the PCF profile provides an unambiguous definition of the structure to be calculated. It is, of course, an idealization of the fiber cross-section as obtained from its SEM image shown in Fig. 1. This is a natural representation for the FEM and the SMT, which make use of the analytically defined locations of glass-air boundaries. On the other hand, in the PWM and the LFM the structure is represented by a superposition of continuous basis functions, which although more adequate for the description of general, graded index profiles, is less suitable for handling abruptly changing profiles. As it was shown by the simple analysis presented in section 2.2, using the plane-wave representation, that is still more accurate than the SEM image, introduces the effective index error of the order of 10^-3. Data from Table 1 indicate that the PWM and the LFM indeed converge to different results than the FEM and the SMT, with the differences for the fundamental modes equal respectively to 4×10^-4 and 3×10^-4, which are lower than the estimated error of the PWM. Obviously, this order of inaccuracy needs to be accepted if the calculated results are matched to experiment. On the other hand, this error has a systematic character and is partially eliminated in the birefringence due to differential nature of this parameter. Therefore, in this work we present independent comparisons for the effective indices and differential parameters such as phase and group birefringences, and GVD.

Out of the rigorous methods, the smallest difference is observed between the SMT and the FEM, while the largest between the SMT and the LFM, systematically for all the simulated parameters. As one should expect, the discrepancies between the different methods increase against wavelength. This effect is intuitively well understood, because for longer wavelength the modes spread deeper into the glass bridges supporting the core and, as a consequence, the calculation results become more sensitive to possible differences in the representation of the input fiber geometry used in respective methods.

 

Fig. 6. Spectral dependence of phase effective index (a) and group effective index (b) of x-polarized LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ spatial modes calculated using FEM.

Download Full Size | PDF

Table 2. Results on the PCF characteristics obtained with the different methods at λ=0.4 µm (first row) and 1.55 µm (second row).

View Table | View all tables in this article

 

Fig. 7. Group velocity dispersion calculated for x-and y-polarized LP ₀₁ (a–b), ${\mathit{\text{LP}}}_{\mathit{11}}^e$ (c–d), and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ (e–f) modes using FEM, SMT, PWM, LFM, and equivalent fiber method EFM.

Download Full Size | PDF

In case of n, the differences between the SMT and the FEM at λ=1.55 µm are equal to 1.1×10^-5, 1.6×10^-5, and 5×10^-5 for the LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes, respectively. Out of the rigorous methods, the differences between the SMT and the LFM are the largest and at λ=1.55 µm are equal to 1.6×10^-3, 4×10^-3, and 4×10^-3 for the LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes, respectively. The discrepancies between the SMT and the PWM for the phase effective index n are systematically about 4 times lower than the ones between the SMT and the LFM. For the group effective index n _g , the dissimilarities are of the same order (SMT-FEM: 7×10^-6, 6×10^-6, and 4×10^-5 for the LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes, respectively; SMT-LFM: 1.1×10^-3, 3×10^-3, and 3×10^-3 for the LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes, respectively). It is interesting to note, that the results obtained with the approximate EFM are in reasonably good agreement with the rigorous methods. In Fig. 7, we show the results of the calculations of the GVD for all the spatial and polarization modes. Similarly as for the phase and the group effective indices, the best consistency is observed between the SMT and the FEM, while the largest difference is between the SMT and the EFM. Nevertheless, the agreement in the zero group velocity wavelength, is very good including the approximate EFM. In particular, for the fundamental doublet, the GVD passes through zero at approximately 914 nm (${\mathit{\text{LP}}}_{\mathit{01}}^x$ ) and 907 nm (${\mathit{\text{LP}}}_{\mathit{01}}^y$ ) according to the FEM simulations, while the EFM gives 897 nm and 884 nm, respectively. Thus, the discrepancy in the zero group velocity wavelength for the fundamental mode does not exceed 2%. The largest difference is found between the FEM and the EFM reaching 2.8% for ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes.

 

Fig. 8. Phase and group modal birefringence calculated for LP ₀₁ (a), ${\mathit{\text{LP}}}_{\mathit{11}}^e$ (b), and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ (c) modes using FEM, SMT, PWM, LFM, and equivalent fiber method EFM.

Download Full Size | PDF

The form anisotropy of the considered PCF gives rise to a noticeable birefringence. In Fig. 8, we present the results of the simulations for the phase and the group modal birefringence defined respectively as:

and

where β _x and β _y are the propagation constants of the guided modes whose degeneracy is removed by the lack of hexagonal symmetry of the PCF. The birefringence simulations are carried out for the LP ₀₁ , ${\mathit{\text{LP}}}_{\mathit{11}}^e$ , and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ spatial modes. Because modal birefringence is a differential parameter, we observe much larger discrepancies between the results obtained with the different methods. Similarly as for the other PCF characteristics, the smallest difference is observed between the SMT and the FEM, while the largest difference is between the EFM and the LFM. The general birefringence characteristics such as the dependence of Δn upon mode order, the negative sign of Δn for the ${\mathit{\text{LP}}}_{\mathit{11}}^e$ mode, and the chromatic dispersion of Δn obtained with rigorous methods are in qualitative agreement with the predictions of the EFM. Nevertheless, the relatively large differences between the EFM and the other methods prove that the glass bridges supporting the PCF core contribute significantly to the overall fiber birefringence. This effect is especially well pronounced for the ${\mathit{\text{LP}}}_{\mathit{11}}^o$ spatial mode, which has the highest overlap coefficient between the field profile and the glass bridges supporting the core, see Fig. 5. The lowest discrepancies between Δn and Δn _g calculated using the rigorous methods and the EFM and therefore the lowest impact of the glass bridges, is observed for the ${\mathit{\text{LP}}}_{\mathit{11}}^e$ mode. This effect is most probably related to asymmetrical shapes of the glass bridges. As can be seen in Fig. 5b, the glass bridges oriented parallel to x-axis are narrower than the others, which results in a decreased perturbation coefficient for the ${\mathit{\text{LP}}}_{\mathit{11}}^e$ mode.

 

Fig. 9. Confinement losses for x-polarized (a) and y-polarized (b) LP ₀₁ , ${\mathit{\text{LP}}}_{11}^e$, and ${\mathit{\text{LP}}}_{\mathit{11}}^o$ modes calculated as a function of wavelength using FEM with PML boundary conditions (solid lines) and SMT (dashed lines).

Download Full Size | PDF

The confinement losses in the PCF are governed by the coefficients a _j (q) in Eq. (4), which quantify the efficiency of coupling between the guided and the radiation modes of an unperturbed fiber. Based on this argument, it is natural to expect that thin silica bridges supporting the fiber core in our PCF structure can only lead to low confinement losses. Indeed, the simulations carried out using FEM with PML boundary conditions and SMT demonstrate (Fig. 9) that the magnitude of the confinement losses for the doublet of LP ₀₁ modes is less than 10^-6 dB/m within the entire visible and near-infrared range down to 2.3 µm. The confinement losses increase for the higher order PCF modes, which is a direct consequence of Eq. (6) for the coupling coefficients a _j (q), showing that the efficiency of coupling is controlled by the overlap of the spatial field profiles in the unperturbed guided and radiation modes with the contours of the fused silica bridges in the PCF. Our predictions for the confinement losses in the considered PCF agree within an order of magnitude with the results of earlier numerical studies [35–37] of confinement losses in several types of PCFs as a function of the number of air holes surrounding the fiber core and the ratio of the hole diameter to the pitch distance.

5. Conclusions

Four rigorous fully vectorial numerical methods (FEM SMT, PWM, and LFM) and approximate equivalent fiber method (EFM) are used in this work to model and better understand the fine details of the field intensity profiles, the field polarization patterns, the phase and the group effective indices, the birefringence, the dispersion properties, and the confinement losses of several spatial modes supported by a multimode photonic crystal fiber with elliptically deformed core and noncircular air holes with high fill factor. With a spline representation of the fiber cross section used as input data, the four rigorous approaches are shown to converge well in their predictions for the phase and the group effective indices, the birefringence and the dispersion characteristics, especially for the lower order modes. For all the analyzed parameters the smallest differences were observed between the SMT and the FEM, while the largest differences between the SMT and the LFM. As explained in Section 4, the residual differences between these rigorous approaches that in case of effective index of fundamental mode are lower than 4×10^-4, are the result of inaccuracies in the representation of the spline-defined fiber geometry in different basis functions. The splines are natural representation for the FEM and the SMT, which make use of the analytically defined locations of glass-air boundaries. In the PWM and the LFM the fiber structure is represented with continuous functions, which are more adequate for the description of graded index profiles.

The results obtained with the approximate EFM are in relatively good agreement with the rigorous methods. For the phase and the group effective indices, the differences between the SMT and the EFM at λ=1.55 µm are of the order (1÷6)×10^-3 depending on the mode order, and are only about two times greater than the differences between SMT and LFM. The discrepancy in calculating the zero group velocity wavelength is also small and ranges from 2% for the LP ₀₁ up to 2.8% for the ${\mathit{\text{LP}}}_{\mathit{11}}^o$ mode. This proves that approximate equivalent fiber methods can be effectively used for predicting the basic characteristics on nonlinear PCFs. More significant discrepancies are observed for the phase and the group modal birefringences, which are related to the differential character of those parameters. The differences between the values of Δn calculated using the FEM and the EFM, observed even for the fundamental mode LP ₀₁ , prove that the glass bridges supporting the PCF core contribute significantly to the overall fiber birefringence. This effect is especially well pronounced for the ${\mathit{\text{LP}}}_{\mathit{11}}^o$ spatial mode, which has the highest overlap coefficient between the field profile and the glass bridges supporting the core.

We have also applied qualitative arguments to illustrate the physical mechanisms by which one can tailor the field profiles, the dispersion, the birefringence, and the confinement losses of the guided modes. In our qualitative treatment, the refractive index profile of the real PCF is viewed as a perturbed index profile of an equivalent elliptical core fiber with known solutions for the guided and the radiation modes. The dispersion properties, the birefringence, the field profiles, and the confinement losses of the PCFs belonging to this class can be qualitatively understood and sometimes well quantified by analyzing the spatial overlap between the eigenmodes of the undistorted fiber and the contours of the correction to the refractive index transforming the idealized elliptical core fiber into a PCF.