FDTD analysis of modal dispersive properties of nonlinear photonic crystal fibers

This paper presents a full-wave electromagnetic analysis of soft-glass photonic crystal fibers developed for the generation of supercontinuum based on third-order nonlinearity. It is shown that a two-dimensional finite-difference time-domain method for guided problems provides results very similar to the measurement data of real fiber structures, enabling the reduction of costly hardware prototyping, thus, opening the way for the application of FDTD to the modeling of nonlinear optical processes.


2005)
, and finite-difference frequency-domain (FDFD) method, (Zhu and Brown 2002). Since dispersive properties of PCFs, critical for nonlinear phenomena, are very sensitive to geometrical properties of a structure, it is essential to have access to a highly reliable EM simulation method in order to reduce costly hardware prototyping to the necessary minimum.
In this paper, a two-dimensional FDTD method developed for guided problems (Gwarek et al. 1993) is applied in a design cycle of PCFs (Salski et al. 2010). The performance of the method is validated against the FDFD method as another widely applicable tool in computational electromagnetics. Both approaches are compared in terms of computational effort and accuracy of obtained results. Eventually, the results are verified against measurement data obtained for microstructured optical fibers (Joannopoulos 2008).

Electromagnetic analysis of PCF
PCF with a solid glass core and the cladding composed of a hexagonal air holes lattice is considered (see Fig. 1). Assuming that the geometry of the fiber does not change along a fiber's principal axis, the problem can be reduced to a vector two-dimensional one (V2D), solvable at a plane transverse to the propagation with an analytically imposed longitudinal phase shift (Gwarek et al. 1993): where λ 0 stands for the wavelength in free space and n eff is an effective refractive index of the considered mode. Consequently, computational effort of the EM analysis can be substantially reduced. In addition, since the analysis is restricted to the fundamental mode with a priori known symmetries, the FDTD model can be reduced to a quarter of the fiber's cross-section with electric and magnetic symmetry conditions imposed, as indicated in Fig. 1. The model is truncated with a perfectly matched layer (PML) surrounding the fiber (Berenger 1996). Before the analysis is started, material properties need to be determined and properly represented with the models available in FDTD. The fibers investigated in this paper are made of Schott glasses, the refractive index of which can be represented by the Sellmeier equation: where n denotes the refractive index of the material, λ is wavelength in microns, B 1,2,3 and C 1,2,3 are experimentally determined Sellmeier coefficients.
Since the Sellmeier equation does not have its direct implementation in the FDTD method, a triple-pole Lorentz model is applied instead: where ε ∞ is optical relative permittivity, ε S is static relative permittivity, f p1,2,3 stand for poles' frequencies, v c1,2,3 represents relaxation frequencies, and A 1,2,3 denote weight coefficients of the corresponding dispersive poles.
As it can be noticed, relaxation frequencies v c1,2,3 occurring in the Lorentz model contribute to a non-zero imaginary part of the permittivity model as in (3). It seems to be in contrast with (2), which is purely real. However, according to the Kramers-Kronig theorem (Landau and Lifshitz 1960), frequency dispersion of a real part of a complex function imposes an imaginary part to be non-zero as well. Although imaginary components are neglected in the Sellmeier equations, it is not allowed to do so in FDTD, which explicitly solves Maxwell curl equations in any causal system, such as the one represented with the Lorentz permittivity model as in (3). Otherwise, an FDTD simulation would have become unstable. For those reasons, relaxation frequencies are set non-zero but small enough (v c1,2,3 = 10 −4 f p1,2,3 ) to assure that the loss factor is negligible. Table 1 shows Sellmeier coefficients and the corresponding Lorentz parameters of two glasses, namely TWNN16 and PBS517 (Stępień et al. 2011). Once the material properties are properly represented, total fiber dispersion can be determined. Total fiber dispersion is comprised of material dispersion and waveguide dispersion. For proper modeling of the fiber both parameters have to be taken into account. Several works previously studied this task (Zhu and Brown 2002;Koshiba and Saitoh 2001;Ferrando et al. 1999).
For total dispersion calculations, several FDTD simulations have to be carried out for a varying longitudinal phase shift β f as in (1), imposed in the V2D FDTD algorithm. In each simulation run, the structure shown in Fig. 1 is excited with a point source driven with the Kronecker delta to evenly cover the investigated spectrum. Subsequently, the Fourier transform of a current injected by the source, as shown in Fig. 2a, is computed to detect resonances indicating guided modes. Those frequencies (or wavelengths) allow determining effective permittivity n eff of the modes for the imposed β f imposed in (1). Figure 2b shows an electric field distribution of the fundamental mode occurring at f = 300 THz for β f = 2,1692 rad/m.
Once the spectrum of the effective refractive index n eff is computed, in the subsequent step group velocity dispersion can be computed as: A similar procedure has to be called if the photonic crystal fiber is supposed to be investigated with the FDFD method (Zhu and Brown 2002). A major difference in respect to FDTD is that a single FDFD simulation provides the solution for a single frequency only, so it has to be executed several times, thus, increasing computational effort.

Ideal PCF structure
A series of EM simulations for a various set of fiber parameters and types of glasses have been undertaken with both FDTD and FDFD methods implemented in QuickWave and Mode Solutions commercial electromagnetic solvers, respectively. Figures 3 and 4 show calculated characteristics of effective refractive index and dispersion of the fundamental mode propagating in PCFs, which are made of PBS517 and TWNN16 glasses. In both cases, the lattice constant of a = 1µm is assumed, whereas the filling factor d/a varies from 0.4 to 0.8. Electromagnetic analysis is undertaken on the Intel Core i5-2410M CPU 2.30 GHz platform. Computing time of full range characteristics with FDFD is around 2 h, while in the case of the FDTD model, computation of the whole characteristic takes ca. 1.5 h. It can be noticed that both methods provide very similar results with a regular discrepancy, although FDFD gives slightly lower values of the effective refractive index (see Fig. 3). The discrepancy increases with the filling factor d/a and with the wavelength λ. It can be also noticed that the application of TWNN16 glass results in more flattened dispersion D in the range from λ = 1 µm up to λ = 2 µm. However, changing the geometry of the photonic crystal cladding, in terms of the filling factor d/a, allows wider adjustment of the dispersion amplitude in the aforementioned spectrum range, which is essential to the efficient supercontinuum generation (compare Fig. 3b with Fig. 4b).
In the next Section, real PCF geometry will be considered and the corresponding computational results will be validated against measurement data, thus, enabling the validation of the applied EM modeling methods.

Real PCF structure
Photonic crystal fiber is developed using a standard fiber drawing technique available at the Institute of Electronic Materials Technology in Poland (see Fig. 5a). The fiber is made of highly nonlinear lead-bismuth-galate PBG08 glass dedicated for supercontinuum generation . Dispersion properties of the fiber are measured with the Mach-Zhender interferometer . The SEM image of the nonlinear PCF has been used as the pattern for the generation of the FDTD and FDFD simulation models.
As it can be seen in Fig. 5c, the measured and calculated dispersion characteristics are in decent agreement, indicating that both numerical methods give reasonable results. However, if the zero dispersion wavelength (ZDW) has to be determined, the FDTD method gives a result much closer to the experimental data. In this case, ZDW is 1,356 and 1,492 nm according to results of the FDTD and FDFD method, respectively, while the measurements indicate that ZDW = 1,382 nm. It means that discrepancy is at the level of about 1.9 and 8.0 % for FDTD and FDFD solutions, respectively. Since determination of ZDW is critical for supercontinuum generation in nonlinear PCFs (Buczynski et al. 2011), the use of the FDFD method can bring a meaningful inaccuracy. Discrepancy between experimental and modeled result is mostly related to credibility of SEM image of the fiber. Since fiber is covers with layer of gold/palladium before the SEM image is taken, diameters of air holes might be slightly different on SEM image and in real fiber. Although the difference is very small it can noticeably influence the dispersion curvature, since dispersion is proportional to the second derivative of effective refractive index (see Eq. 2).

Conclusions
Application of the FDTD method to the computation of modal dispersive properties of photonic crystal fibers has been addressed. The accuracy of the proposed method was verified against the FDFD method and against the experimental results. Due to a frequency-domain approach, the extraction of the effective refractive index with FDFD is based on tracing of a selected mode versus wavelength. On the contrary, FDTD allows tracing all the modes simultaneously, which can be advantageous when modes coupling becomes of our interest. In general, both methods give similar results and can be successfully applied to the modeling of photonic crystal fibers.