{'fraction_non_alphanumeric': 0.042586442664582926, 'fraction_numerical': 0.022953701894901347, 'mean_word_length': 4.320083138477527, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of graded-index optical fibers by the spectral parameter power series method

Spectral parameter power series (SPPS) method is a recently introduced technique for solving linear differential equations and related spectral problems. In the present work we develop an approach based on the SPPS for analysis of graded-index optical fibers. The characteristic equation of the eigenvalue problem for calculation of guided modes is obtained in an analytical form in terms of SPPS. Truncation of the series and consideration in this way of the approximate characteristic equation gives us a simple and efficient numerical method for solving the problem. Comparison with the results obtained by other available techniques reveals clear advantages of the SPPS approach, in particular, with regards to accuracy. Based on the solution of the eigenvalue problem, parameters describing the dispersion are analyzed as well.


Introduction
Analysis of graded-index cylindrical waveguides, in particular of optical fibers, presents considerable mathematical and computational difficulties. The main differential equation involved is singular, the fact which restricts application of purely numerical techniques. The typically used approach is based on the asymptotic WKB approximation (Wentzel-Kramers-Brillouin). It can offer analytical relations between the mode-propagation parameters (see, e.g., [21], [23,Sect. 3.7]). However, it is well known (see, e.g., the discussion in [20]) that they are not accurate enough for the practical study of optical fibers.
We develop a completely new approach for analysis of the computationally difficult spectral problems involved. It is based on the spectral parameter power series (SPPS) representation for solutions of second-order linear differential equations of the Sturm-Liouville type. The SPPS representation for regular Sturm-Liouville equations was first obtained in [18] and applied to numerical study of spectral problems in [19]. It was used in a number of applications (we refer to the review [16]). In particular in [7] the SPPS method was applied to the study of wave propagation in layered media. In [8] the SPPS representation was obtained for solutions of the singular differential equations belonging to the class of perturbed Bessel equations. The main differential equation of the graded-index cylindrical waveguides is of that type. In the present paper we use the main result of [8] to develop a numerical method for analysis of graded-index optical fibers. The characteristic equation equivalent to the spectral problem for calculation of guided modes is obtained in an ex-plicit analytical form. The numerical method consists in approximating this equation and finding zeros of the approximate characteristic function.
The paper contains several numerical tests. In the case of the well studied parabolic profile we compare our results with the available exact solution, with the WKB approach as well as with the results from [20] where the finite element method was applied with the results of the WKB approximation used as the initial guess. Our method gives substantially more accurate values.
We show that the SPPS method presented here allows one a quick and accurate analysis of graded-index optical fibers and can be used as an efficient tool for their design.

Mathematical formulation of the problem
The basic wave equation governing the wave propagation in graded-index fibers has the form (see, e.g., [23,Sect. 3.7]) where the knowledge of the solution ψ allows one to compute the corresponding electromagnetic field, n(r) is the radial refractive index profile, k = 2π/λ is the vacuum wave number, β is the propagation constant, and m is a mode parameter given by As in most practical fibers the refractive index varies in the core but is constant in the cladding, it is usually convenient to solve the wave equation in the core and cladding separately and to match those solutions at the core-cladding boundary. If the solutions in the core and cladding are denoted by ψ(r) and ψ clad (r), respectively, the boundary conditions at the interface r = a under the weakly guiding approximation (see, e.g., [23, p. 94]) are given by ψ(a) = ψ clad (a) and dψ dr r=a = dψ clad dr r=a .
We note that for the method presented here it is not essential that the boundary conditions be that simple. For example, they could be dependent on k.
Introducing the function U (r) = √ rψ(r) we write equation (2.1) in the form Finally, introducing the new independent variable x = r/a we obtain an equation for the function The boundary conditions (2.2) take the form The solution for the cladding (x > 1) has the form where K m is the modified Bessel function of the second kind, the constant n 2 is the value of the refractive index in the cladding, and C is an arbitrary constant. The bounded solution u in the core is unique up to a multiplicative constant. Thus, from the boundary conditions (2.4) the following characteristic equation of the problem is obtained which then admits the form Pairs of values of β and k which satisfy this characteristic equation together with the light propagation condition k 2 n 2 2 < β 2 < k 2 n 2 1 := k 2 max 0≤x≤1 n 2 (ax) (2.6) correspond to guided modes in the fiber and are the main object of computation. As we show below, the SPPS approach allows one to approximate directly the characteristic equation and to solve it with a considerable accuracy and speed.
Remark 2.1. Note that the characteristic equation (2.5) remains valid also for the case of absorbing media, i.e., when the refractive index n(r) has non-zero imaginary part. The only change is in the propagation condition (2.6), it should be written as Re n 2 (ax).

Analysis of dispersion
The solution of the characteristic equation (2.5) gives us the dependence of k on β (or, since k = 2π/λ = ω/c, the dependencies of λ and ω on β). A way to employ the obtained information is in the analysis of dispersion. The modal description of dispersion is related to the different mode indices (or group velocities) associated with different modes. The group velocity associated with the fundamental mode is frequency dependent because of chromatic dispersion. As a result, different spectral components of the pulse travel at slightly different group velocities, a phenomenon referred to as group velocity dispersion (GVD) which has two contributions: material dispersion and waveguide dispersion [2]. Group velocity v g is defined as [5] v g = dω/dβ (3.1) where ω stands for the frequency. If ∆ω is the spectral width of the pulse, the extent of pulse broadening for a fiber of length L is governed by [2] ∆T = d dω where ∆T stands for the differential group delay and the parameter β 2 = d 2 β/dω 2 is known as the GVD parameter. In terms of the range of wavelengths ∆λ emitted by the optical source, and by using ω = 2πc/λ and ∆ω = (−2πc/λ 2 )∆λ, (3.2) can be written as where the dispersion parameter D, expressed in units of ps/(km-nm), is given by Interest in studying fibers with different refractive index profiles comes from the fact that tayloring such profiles the dispersion, which is one of the main impairments for optical communications, can be manipulated. Fibers where the dispersion parameter can be reduced at a certain wavelenght, flattened for a wide range of wavelenghts or even made highly negative are used for aplications concerning long haul communications, wavelenght division multiplexing and dispersion compensation, among others (see, e.g., [3], [22,Chapter 4]).

Preliminary facts and SPPS representations for perturbed Bessel equations
The main equation where m ∈ N, q and r are known functions (which in general can be complex valued), µ is a spectral parameter. The equation belongs to the class of perturbed Bessel equations and was studied in a considerable number of publications (e.g., [6], [8], [9], [13], [17], [24]). In what follows we always assume that q and r are piecewise continuous functions on the whole segment of interest with at most a finite number of step discontinuities. The SPPS method applied to (4.1) consists of two steps. On the first step one needs to construct a nonvanishing on (0, a] bounded solution u 0 of the equation In [8] it was proved that if q(x) ≥ 0, x ∈ (0, a] then such solution u 0 exists and has the form where the functions Y (j) are defined recursively as follows  The series (4.3) converges uniformly on [0, a]. When q does not satisfy the nonnegativity condition the series (4.3) still defines a bounded solution of (4.2) but in general not necessarily nonvanishing on (0, a]. The derivative of u 0 has the form On the second step the solution u 0 is used for constructing the unique (up to a multiplicative constant) bounded solution of (4.1) for any value of µ. This solution has the form where the functions X (j) are defined recursively as follows The series (4.6) converges uniformly on [0, a]. The first derivative of u is given by where the series converges uniformly on an arbitrary compact K ⊂ (0, a]. In (4.6) the solution is represented in the form of a power series with respect to the spectral parameter µ. Such representations were called in [19] the spectral parameter power series (SPPS). The representation (4.6) is based on a particular solution of equation (4.1) for µ = 0. Similarly to the Taylor series, the approximations given by the truncation of the series (4.6) and (4.8) are more accurate near the origin, while the accuracy deteriorates with the increase of the parameter µ, see [8,Example 7.7]. In [8,19] it was mentioned that it is also possible to construct the SPPS representation of a general solution starting from a non-vanishing particular solution for some µ = µ 0 . Such a procedure is called a spectral shift and consists in the following. Equation (4.1) can be written in the form where µ 0 is an arbitrary number. This equation is again of the form (4.1) but with the spectral parameter Λ := µ − µ 0 . The SPPS technique applied to this equation leads to a representation of the solution in the form of a power series in terms of Λ and allows one to improve the accuracy of the approximations given by the truncated series (4.6) and (4.8) near the point µ 0 . The spectral shift has already proven its usefulness for numerical applications [19], [15], [8] and is used below in numerical computations. The sufficient condition q(x) − µ 0 ≥ 0 mentioned above imposes an upper bound for the possible choices of the parameter µ 0 . Allowing µ 0 to be complex valued can solve the problem. Assume that q is real valued and r(x) > 0 almost everywhere on (0, a] (which is sufficient for the scope of the paper). Then every particular bounded solution of equation (4.1) having µ such that Im µ = 0 is non-vanishing on (0, a]. Indeed, the equality u(x 0 ) = 0, x 0 ∈ (0, a] contradicts the well-known property (see, e.g., [26,Chapter 10.4]) of the operator τ u = −u ′′ + m 2 −1/4 x 2 + q(x) u with a suitable domain to be self-adjoint in the Hilbert space L 2 ((0, x 0 ], r) and have only real eigenvalues.

The SPPS form of the characteristic equation
Let the parameter β be fixed. Writing equation (2.3) as we obtain an equation in the form considered in Section 4 with the spectral parameter µ = a 2 k 2 and having known particular non-vanishing solution u 0 = c √ xI m (aβx). However the direct application of the SPPS method for the obtained equation leads to certain difficulties. First, for the practically used fibers the value aβ is quite large (in the examples below it is within the range 20−300) leading to a rapidly growing particular solution u 0 . Second, since we are looking for the values of the parameter k satisfying (2.6), we cannot take advantage of the property of the SPPS representation to be the most accurate near the origin. To overcome these difficulties we apply the spectral shift technique. Let us assume that the refractive index n is real valued. By n 1 we denote the maximum of n in the core. Note that n 1 > n 2 . The main equation (2.3) can be written in the following form where k 0 = β/n 1 . Then q(x) := a 2 β 2 − k 2 0 n 2 (ax) ≥ 0, x ∈ (0, 1]. Introducing the notations µ := a 2 k 2 − k 2 0 and r(x) := n 2 (ax) we find that equation (5.1) is an equation of the form (4.1) with q satisfying the nonnegativity condition from the preceding section. Hence the solution u 0 constructed as explained above does not have other zeros on [0, 1] except at x = 0 and the bounded solution of (5.1) has the form (4.6).

Numerical implementation and examples
Based on the results of the previous sections we can formulate a numerical method for computing the guided modes of the fiber. We suppose that some range for the parameter β is given, say β ∈ [β 1 , β 2 ]. If instead a certain segment of wavelengths [λ 1 , λ 2 ] is given, one can find the range of β's from (2.6):

Find its zeros satisfying propagation condition
6. If the interval for the spectral parameter µ in (5.1) defined by (6.2) is large (e.g., max µ > 400), perform several steps of the spectral shift technique to improve the precision. A non-vanishing particular solution on each step can be obtained using Remark 4.1.

7.
Combine the values of parameter k found for different values of β to obtain dependencies k(β) for different propagation modes.
Before considering numerical examples let us explain how the numerical implementation of the SPPS method was realized in this work. All calculations were performed with the help of Matlab 2013b in the double precision machine arithmetics in a PC with Intel i7-3770 processor. The formal powers X (j) were calculated using the Newton-Cottes 6 point integration formula of 7 th order, see, e.g., [10], modified to perform indefinite integration. We choose M equally spaced points covering the segment of interest and apply the integration formula to overlapping groups of six points. It is worth mentioning that for large values of the parameter m a special care should be taken near the point 0, because even small errors in the values of X (2j−1) after the division by u 2 0 ∼ x 2m+1 lead to large errors in the computation of X (2j) on the whole interval [0, 1]. To overcome this difficulty we change the values of X (2j−1) in several points near zero to their asymptotic values. This strategy leads to a good accuracy. The computation of the first 100 to 200 formal powers proved to be a completely feasible task, and even for M being as large as several millions the computation time of the whole set of formal powers is within seconds. In the presented numerical results we specify, among others, two parameters: N is the number of computed formal powers X (j) , and M is the number of points taken on the considered segment for the calculation of integrals. Using the calculated formal powers we evaluated the approximation of the characteristic equation (5.2), constructed a spline passing throw the obtained values and found its zeros using the Matlab command fnzeros.

Mode characterization of a graded index fiber
As a first case of study, we present the computation of the propagation constants β for a graded index multimode fiber with a refractive index given by where the relative refractive-index difference between n 1 and n 2 is defined as For α = 2 the exact solution of the problem is known allowing us to verify the accuracy of the proposed method. The specific values that were used in order to be able to compare our results with the ones in the literature and with an exact solution were n 1 = 1.462, n 2 = 1.447, a = 25 µm, and λ = 0.78 µm. The parameters associated to the numerical implementation of the SPPS method included the approximation of the functions with M = 100001 points and using N = 100 formal powers X (j) for the truncated series. For this example it was natural to apply Remark 5.1. Since the range of the new spectral parameter µ = a 2 (β 2 0 − β 2 ) defined by (2.6) is approximately [0, 1700], i.e., quite large, we applied 8 steps of the spectral shift. As a result, on each step the range of the spectral parameter was less than 200. The non-vanishing solutions were obtained accordingly to Remark 4.1, we used spectral shifts µ n = (200 + 4i)n. On each step a 1000 points mesh was used to approximate the characteristic equation.
Our program has found 121 propagation constants ending with m = 20, the overall computation time was about 3 minutes. The number of modes is in a good agreement with that obtained from the WKB approximation [23, 3.7.2]. Indeed, taking into account that there are 2 conventional modes corresponding to m = 0 and 4 corresponding to m ≥ 1 (see [23, 3.4 Table 1 (c) we mention that the propagation constants obtained for m = 0 were the most inaccurate. Such phenomenon has already been observed in [8,Example 7.5] and the accuracy can be improved by taking more points to represent the formal powers X (j) . The accuracy of the obtained propagation constants for larger values of m was excellent even taking less points representing the formal powers. It is also remarkable that the propagation constants conform clusters so that the modes in the same group will propagate with very similar constants. We illustrate this additionally on Figure 1 where the graphs of the absolute value of the characteristic function (2.5) are presented for different values of m.

Dispersion curves of LP modes
The SPPS method allows us not only to find the propagation constants for a particular wavelength but also to analyze the behavior of the propagation constants for a range of wavelengths. For a numerical illustration we considered the profile (6.3) for α = 1 (the so-called triangular profile) with n 1 = 1.462, n 2 = 1.447, a = 12.5 µm. On Figure 2 we present the graphs of the normalized propagation constant versus normalized frequency V = ak n 2 1 − n 2 2 for different propagation modes. The parameters of the SPPS implementation were M = 10001, N = 100 and we used 140 different values of V to calculate the propagation constants. We refer the reader to [12] where the normalized cut-off frequency V ≈ 4.381 of the single-mode operation was computed for the triangular profile fiber. Our results agree with this value.

The group velocity and the waveguide dispersion
The SPPS method was used to evaluate the group velocity (3.1) and the dispersion coefficient (3.3) for the triangular profile (with the refractive index defined by (6.3) with α = 1). The parameters of the fiber were taken from [3] and were the following: n 1 = 1.527, n 2 = 1.5094395, a = 3.2 µm. For the computation we used 1.00 ≤ λ ≤ 1.60 µm, 1000 values for k, 100 values for β, M = 10001 and N = 100. The elapsed time was 0.56 minutes. The group velocity is shown on Figure 3, and the dispersion parameter is shown on Figure 4. In both graphs the refractive index profile of the fiber was included in a subplot.

Simultaneous analysis of material and waveguide dispersion
The refractive index depends on the wavelength. Such dependence is usually described by the Sellmeier equation, see, e.g., [1, (7.126)] or by its generalizations allowing graded-index fibers where dopant concentrations depend on the radius, see [11], [14], [25]. The SPPS method applied as described in Remark 5.1 can be used to study the combined material and waveguide dispersion (chromatic dispersion).
As an example we considered a dispersion-flattened triple clad fiber (fiber 2 from [4]). The dopant concentrations of the fiber are the following: core -9.1 m/o P 2 O 5 , cladding 1 -13.5   Dispersion coefficient is shown on Figure 5 , and the group velocity is shown on Figure 6. In both graphs the refractive index profile of the fiber was included in a subplot. As can be seen from the graph, the dispersion is within ±1 ps/(km·nm) for 1.382 ≤ λ ≤ 1.586 µm.

Conclusions
The SPPS method is developed and applied to problems of wave propagation in graded-index optical fibers. The numerical examples presented show that the method provides a remarkable accuracy achievable within seconds. It shares with the widely used WKB approach the possibility to work with an analytic representation for the solution and for the characteristic function of the problem, offering at the same time the accuracy superior to other available purely numerical methods such as finite element method. We believe that the SPPS approach will become a standard tool for analysis and design of optical fibers and other inhomogeneous cylindrical waveguides.