1. Introduction

As part of nonlinear optical effects, the optical parametric processes are playing an important role in the proposal of a new generation of optical devices. They are called parametric because the processes are originated from light-induced modulation of a medium parameter such as refractive index [1] or, equivalently, the propagation constant β. As this is frequency dependent, mathematically the effects of medium dispersion can be taken into account by expanding β of the considered mode in a Taylor series around a central frequency${\omega }_0$:

From Eq. (1), the behavior of β in the neighborhood of ${\omega }_0$ depends on the ${\beta }^{(i)}$’s at ${\omega }_0$, which are defined in Eq. (2) and are the so-called Higher-Order Dispersion Parameters (HODP).

The parametric processes’ efficiency in a given device depends on phase matching, which in turn depends on HODP’s values. Therefore, a reliable method to compute these parameters would be of great value to the optical waveguide designer. Furthermore, from a designing view point, HODP up to sixth-order may be considered adequate to model satisfactorily the nowadays optical devices and, probably, those that will emerge in a near future, therefore, those HODP have been the focus of the present analysis.

Given the importance that the knowledge of ${\beta }^{(i)}$’s represents to optical parametric processes, studies on the numerical precision required to obtain these parameters with an adequate degree of reliability seem to be absent. Thereby, this work constitutes, to the best of our knowledge, the first systematic study of HODP’s numerical calculation errors reported in the literature.

In order to guarantee that a numerical strategy can actually get adequate results, it must undergo a rigorous numerical assessment process through the analysis of several key geometries. To the study presented here, the numerical strategy must be able to provide accurate enough${\beta }^{(0)}$’s values (as output of one of its intermediate stages) to ensure that, even after they undergo through a differentiation process, HODP values thus obtained are sufficiently accurate to permit to build up optical devices based on its theoretical prediction. In fact, the knowledge of the code’s accuracy draws a line that defines its limitations and applicability.

Therefore, this study will show, through a numerical assessment process, that the numerical strategy briefly described in Section 2 is capable of providing the required numerical accuracy in general optical waveguides’ HODP calculations. Moreover, although the main building blocks of the strategy are methods well-known in the literature, the conjunction of them itself constitutes a reliable and efficient method to get adequate HODP predictions.

Next, a case study involving Photonic Crystal Fibers (PCF) HODP optimization is presented. PCF are much more complex geometries than the waveguides analyzed during assessment. In order to accomplish this task a Genetic Algorithm (GA) was implemented and coupled to the original numerical strategy. In this case, any proof of the results’ accuracy will not be provided, given the lack of theoretical or experimental results in the literature that provide sufficient accuracy to be used as reference.

According to above discussion, this study is organized as follows: in Section 2 is introduced the numerical strategy; in Section 3 are pointed out some numerical considerations concerning the approach adopted in this work; in Section 4 is performed the assessment of the numerical strategy; in Section 5 is presented a case study of a PCF obtained via GA optimization; in Section 6 the conclusions and final comments are given.

2. Numerical strategy for HODP calculation

The numerical strategy adopted here is constituted by two primarily building blocks: a modal solver and a proper Finite Difference Differentiation Algorithm (FDDA).

The first one is a 2D Full-Vectorial Modal Solver, based on the Finite Element Method (FEM), which is an implementation of the concept presented in [2] with hybrid edge/nodal Finite Elements (FE) [3]. The computational domain’s truncation was made with Perfectly Matched Layers (PML) [4]. In the case of dispersive media the material dispersion is always taken into account in eigenvalue’s calculation.

Additionally, the solver allows a user to make some choices. Rectilinear or curvilinear hybrid FE can be selected and in both cases two subtypes of FE can be chosen: either Constant Tangential/Linear Normal (CT/LN) or Nedelec’s space-conforming Linear Tangential/Quadratic Normal (LT/QN) [3].

A diagram illustrating the numerical strategy’s structure is shown in Fig. 1 .

 

Fig. 1 Diagram of the developed numerical strategy.

Download Full Size | PDF

According to Fig. 1, the FEM Solver supplies the ${\beta }^{(0)}$’s values (eigenvalues) to the FDDA module and this, in turn, is responsible for obtaining the HODP corresponding to each of the waveguide modes. Its current implementation is a result of various numerical considerations, described in Section 3, to ensure numerical accuracy and efficiency.

3. HODP calculation: numerical considerations

3.1 Finite element choice and numerical integration

The LT/QN FE were adopted here due to its superiority over CT/LN ones for calculating ${\beta }^{(0)}$ related to optical waveguides, as stated in [3]. Notice that no study to date has shown that results obtained by LT/QN FE is sufficiently accurate to permit to get a good precision in HODP calculations even after a subsequent application of some differentiation method.

Additionally, when using curvilinear FE the elementary integrals must be calculated numerically. Gaussian Quadrature was employed for this. Tests considering 7 and 13 sampling points per triangle were performed and results showed no significant error difference between these two cases. Thereby, 7 points were adopted in all this study.

3.2 Differentiation method choice

Regarding differentiation, two ways are possible: straight differentiation through a numerical method or work with polynomials making a fitting and then applying analytical differentiation. For experimental data differentiation, both paths would be tortuous due to uncertainties and noise. The second path is usually preferred because numerical differentiation becomes extremely ill-conditioned in this case. However, in the theoretical case (data from simulations) continuous and smooth curves can be obtained, provided that certain precautions are taken. Thereby, the choice of which of those two methods is the best to be applied here would be arguable and beyond this work’s scope.

In this study, finite difference numerical differentiation was chosen because it is an extremely well established method and their formulas are extremely simple to implement computationally. All formulas based on finite differences used in this work were obtained from [5]. Approximation errors resultants of FDDA application are the truncation, rounding and cancellation ones. All of them directly influence the accuracy of results and are intrinsically connected to the differentiation step h. The choice of h is a compromise. Unfortunately, there is no way to estimate in advance what is the best value to be chosen. Here h’s choice was made by several numerical simulations which allowed finding a more appropriate h. Thus, for all simulations was taken$h\cong 8.58\times {10}^{12}rad/s$.

Additionally, the process of obtaining HODP consists, in a global sense, of two main steps: calculation of ${\beta }^{(0)}$ by the modal solver and its subsequent numerical differentiation process. The best opportunity to reduce the error appears in the second step, in which one the number of sampling points from the function whose derivative is wanted can be increased, what reduces the truncation error in principle. However, the overall error does not only depend on the truncation error. In fact, a new sampling point addition to the differentiation formula also implies an accumulation of new characteristic rounding error corresponding to that point. The result is an error accumulation so that, as points’ number is increased, these errors compensate the original gain in accuracy by reducing the truncation error. Thus, choosing the appropriate points’ number to be considered in the formula is also a compromise. Here up to sixth-order derivatives will be calculated, what requires at least seven points’ formulas. Accordingly, all assessment analysis items were performed considering formulas with 7, 9 and 11 points. This procedure has permitted to choose which points’ number should be used in the evaluation of more complex waveguide geometries.

3.3 Considerations about chain rule use

With respect to the differentiation parameter, HODP calculation can be done in two different ways: either straightforwardly use${\beta }^{(0)}$’s numerical values versus ω and obtain the other parameters by derivatives with respect to ω (referenced here as ω-approach), or apply the chain rule in order to write the derivatives with respect to ω as derivatives with respect to λ and use these derivatives’ values to calculate HODP (referenced here as λ-approach). From a mathematical view point, both approaches are identical and produce the same result. However, from a numerical view point, sometimes using the second one can produce less accurate results. This is because, using the ω-approach, we may write,

The assessment tests showed that in most of the cases a sufficient error accumulation to cause significant changes in results’ accuracy was not observed. However, in some cases an accuracy reduction of results was apparent for the higher-order derivatives. Thus, as a precaution, whenever possible, derivatives should be performed with respect to ω.

3.4 Error calculation

Aiming at quantifying how good is the approximation provided by the numerical strategy, a quantitative analysis of the error was made. A quantification that could be applied to all dispersion parameters regardless of its magnitude was chosen: the decimal logarithm of the relative error. Here, it will be designated by the acronym RE (initials from Relative Error) and defined by the equation

4. Assessment

4.1 Rectangular waveguide analysis

In order to make a rigorous numerical accuracy analysis of the adopted numerical strategy, the homogeneous rectangular waveguide analytical case, with vacuum as dielectric, was chosen as reference. Waveguide’s walls were considered Perfect Electric Conductors (PEC). Although this is an idealized situation (no metal optical properties have been considered), it represents an excellent starting point for the analysis presented here. In fact, this structure provides explicit analytical expressions for all${\beta }^{(0)}$’s derivatives with respect to ω and, because of that, this case has been chosen as reference for the assessment process.

A schematic diagram and the results from the approximation error analysis are presented in Fig. 2 . In all simulations the parameters a = 1 μm e b = a/2 were assumed [see Fig. 2(a)]. Due to its symmetry only half waveguide was simulated utilizing rectilinear FE.

 

Fig. 2 Analysis of the rectangular waveguide. (a) Geometry of the guide; (b) RE versus Degrees of Freedom for 1299.41 nm; (c) RE versus Degrees of Freedom for 1702.33 nm; (d) RE versus wavelength for 75191 DOF.

Download Full Size | PDF

After conducting numerous tests no significative accuracy gain was noted when using 9 or 11 sampling points instead of 7 ones in derivatives. Therefore, in the rectangular waveguide’s case was found that 7 points are enough to model with necessary accuracy its HODP. So, the following results were obtained for this case.

HODP’s study was performed in 1.3 – 1.7 μm band (optical communications’ band). Two points close to its extremes were picked up: 1299.41 nm and 1702.33 nm. For these two wavelengths the RE, given in Eq. (5), was evaluated and plotted versus the number of Degrees Of Freedom (DOF) in the mesh. The corresponding results are shown in Figs. 2(b) and 2(c), for 1299.41 nm and 1702.33 nm, respectively. In Fig. 2(d) the RE versus wavelength is presented for 75191 DOF.

At Fig. 2 the relative error tends to increase when increasing derivatives’ orders to be calculated. Accordingly, lower accuracy is obtained for the highest-order derivatives. In general, to build devices with controlled fourth-order dispersion is currently complicated because process’ inherent fabrication errors. However, in some cases in which geometry isn’t very complex, to produce such devices could be possible [6]. So, if a sixth-order controlled dispersion device could be accurately predicted would be a great achievement, because this one will probably constitute a result of a near future optical waveguides’ state of the art manufacturing processes result.

Tables 1 and 2 show, respectively for 1299.41 nm and 1702.33 nm, the RE for 75191 DOF. The largest error is obtained for${\beta }^{(6)}$, which is about 0.23% for 1299.41 nm and about 2.5% for 1702.33 nm. These precision values are considered very good even when considering up to sixth-order dispersion. If only up to fourth-order dispersion is considered, the errors are 0.00% and 0.01% for 1299.41 nm and 1702.33 nm, respectively. Finally, from Figs. 2(b) and 2(c) may be noted that for 10000 DOF a good approximation has already been reached.

Table 1. HOD parameters of the rectangular waveguide at λ = 1299.41 nm. They were used 7 sample points in the derivatives calculation and 75191 DOF in the mesh

View Table | View all tables in this article

Table 2. HOD parameters of the rectangular waveguide at λ = 1702.33 nm. They were used 7 sample points in the derivatives calculation and 75191 DOF in the mesh

View Table | View all tables in this article

4.2 Step-index profile fiber analysis

Carrying on with the assessment process, a structure with curved interfaces and open boundaries is evaluated. For this, the fiber with step index profile proposed in [6] and called F3 was picked up. That fiber is included in the W-profile category. In this case, the semi-analytical model’s results from [6] are taken as reference. A schematic diagram and the results from the approximation error analysis are presented in Fig. 3 .

 

Fig. 3 Analysis of the Fiber F3. (a) Geometry of the W-profile fiber analyzed in [6]; (b) RE versus Degrees of Freedom for 1299.41 nm; (c) RE versus Degrees of Freedom for 1702.33 nm; (d) RE versus wavelength for 75099 DOF.

Download Full Size | PDF

F3′s profile is showed in Fig. 3(a). For this waveguide is not possible to obtain analytical expressions for HODP. However, the solution accuracy of the corresponding eigenvalue’s equation is very high (about machine ε), allowing for results to be used as reference in error calculation. Because numerical data is in tabular form in this case, they must undergo a numerical differentiation process in order to obtain the remaining HODP. This implies that the reference has inherent differentiation process’ errors.

Concerning the number of sampling points in the differentiation formula, the same as observed in the rectangular waveguide case is kept for the fiber with step-index profile under analysis: 7 sampling points are sufficient to obtain the necessary accuracy. Thus, all results from this subsection were obtained for this case.

Due to F3′s symmetry, only one quarter of the fiber was simulated using curvilinear FE. The simulations are also made here in the optical communications’ band. The RE versus DOF for 1299.41 nm and 1702.33 nm is shown in Figs. 3(b) and 3(c), respectively. In Fig. 3(d) the RE versus wavelength for 75099 DOF is displayed. Here the lowest accuracy is also obtained for the highest-order derivatives.

Tables 3 and 4 show, respectively for 1299.41 nm and 1702.33 nm, all numerical values of interest. The maximum error occurs for${\beta }^{(6)}$: about 0.28% for 1299.41 nm and approximately 0.89% for 1702.33 nm. These precision values are considered excellent for the case of sixth-order dispersion. If only up to fourth-order dispersion is considered, the maximum errors are 0.00% and 0.01%, for 1299.41 nm and 1702.33 nm, respectively. Furthermore, analyzing Figs. 3(b) and 3(c), an excellent approximation could have already been obtained for a little more than 20000 DOF. In a global sense, an impressively good approximation was obtained.

Table 3. HOD parameters of the step index fiber at λ = 1299.41 nm. They were used 7 sample points in the derivatives calculation and 75099 DOF in the mesh

View Table | View all tables in this article

Table 4. HOD parameters of the step index fiber at λ = 1702.33 nm. They were used 7 sample points in the derivatives calculation and 75099 DOF in the mesh

View Table | View all tables in this article

5. An inverse problem in PCF

Once evaluated the numerical strategy and shown that it provides reliable values for calculating HODP of conventional waveguides, now more complex waveguides such as PCF can be considered. The high number of degrees of freedom that this guide may display in its structure (materials, hole’s diameter, pitches, etc.) permits great flexibility in obtaining the most diverse requirements for HODP, however, it makes the problem of synthesis quite complex [7]. To address such inverse problem, the use of search and/or optimization methods is essential. Among a variety of methods that can be used for such purposes, the GA has been the most widely used for solving similar problems in PCF.

As a consequence, a conventional GA based on that presented in [8] has been implemented and coupled to the original computational code. In Fig. 4 , there is a flow chart describing the numerical strategy automation to deal with inverse problems. In the GA a user may define constraints to PCF structure (such as materials, number of rings, maximum and minimum dimensions, among others) from which are built up the individuals’ chromosomal representations. After that, the structure given by the chromosome is automatically generated, meshed and used as the input for the modal solver. The output data of this one suffers a post-processing in the FDDA module, from which are obtained the HODP, effective areas, nonlinear coefficients, etc. These outputs are combined into a mono or multi-objective function that constitutes the problem fitness criterion.

 

Fig. 4 Solver automation for the treatment of inverse problems.

Download Full Size | PDF

In order to demonstrate the effectiveness of the scheme presented in Fig. 4, among many possibilities of problems to be solved, a simple case study based on design requirements involving HODP is proposed. Therefore, a PCF configuration that should satisfy design requirements involving D e ${\beta }^{(4)}$ was aimed, where D is given by

A schematic diagram for the initial PCF and the HODP results from the final PCF are presented in Fig. 5 . The starting point is a geometric structure with hexagonal arrangement of air holes as shown in Fig. 5(a). The number of periods was taken as constant and equal to 4, number that was considered suitable to supply sufficient degrees of freedom to make the results reach the goals.

 

Fig. 5 The inverse problem solved. (a) PCF obtained by GA; (b) β⁽¹⁾ and β⁽²⁾ versus wavelength; (c) β⁽³⁾ and β⁽⁴⁾ versus wavelength; (d) β⁽⁵⁾ and β⁽⁶⁾ versus wavelength.

Download Full Size | PDF

Equation (7) represents the Fitness Function (F) adopted,

Additionally, a nonlinear coefficient $\gamma \ge 10W^{-1}/km$ is also desirable for a PCF in this type of application and, in spite of not considering this parameter in the function F, a PCF with SF6 glass from Schott Glass [10] was designed. This type of glass has a nonlinear index (n ₂) higher than pure silica (the n ₂ used was obtained via the formula reported in [11] through the data provided by the manufacturer), which facilitates a higher γ attainment. Additionally, the non-inclusion of γ in F allows reducing the problem’s optimization complexity.

About the number of sampling points used in the FDDA, as in all assessment cases 7 sampling points were sufficient, this number was adopted as a standard for this simulation.

In the simulation carried out to solve this inverse problem was used a PC platform with a 2.4 GHz Intel® Core2 Quad processor and 8GB RAM. The GA started with a population of six individuals. The selection method used was Roulette Wheel and the mutation operator performs a random selection of the ring to be muted when it is applied.

After 312 generations (with a total time of 140 hours considering LT/QN FE), the best individual that was obtained had the following geometric parameters: pitch Λ = 5.07μm and hole diameters d₁ = 0.36 μm, d₂ = 0.64 μm, d₃ = 1.44 μm and d₄ = 1.44 μm. Here it is important to point out that all geometric dimensions were swept with a two decimal places’ discretization (a minimum variation of 0.01 μm in all the dimensions was allowed). The obtained nonlinear coefficient was $\gamma \cong 13W^{-1}/km$ at 1.55 μm and the HODP versus wavelength plots can be seen in Figs. 5(b), 5(c) and 5(d). The results are in good accordance with the established goals, as can be seen from those figures.

6. Conclusions

We presented an efficient numerical strategy that permits an accurate analysis of Optical Waveguides’ HODP. A systematic study of the numerical approximation error that results from HODP calculation was performed for the first time to our knowledge. The results have shown that the strategy proposed here can model accurately HODP of that guides. More specifically, from the assessment process, adopting a $RE\le -10$ for ${\beta }^{(0)}$ and $\Delta \omega \cong 8.58\times {10}^{12}rad/s$, then a finite difference formula of 7 sampling points is enough to ensure a $RE\le -1$ for ${\beta }^{(6)}$ and a $RE\le -3$ for ${\beta }^{(4)}$ in all optical communication’s band.

Additionally, a GA was implemented and coupled to the original computational code to solve inverse problems related to HODP in PCF, completing the numerical strategy by adding the synthesis capability to it. Its efficiency was demonstrated through a simple case study where the results obtained were close to the requested goals.