Reconstruction of Dielectric Constants of Core and Cladding of Optical Fibers Using Propagation Constants Measurements

We present new numerical methods for the solution of inverse spectral problem to determine the dielectric constants of core and cladding in optical fibers. These methods use measurements of propagation constants. Our algorithms are based on approximate solution of a nonlinear nonselfadjoint eigenvalue problem for a system of weakly singular integral equations. We study three inverse problems and prove that they are well posed. Our numerical results indicate good accuracy of new algorithms.


Introduction
In this paper, we investigate the problem of determination of dielectric constants of core and cladding in optical fibers using measurements of propagation constant.Such kinds of problems arise in determination of electromagnetic parameters in nanocomposite materials in nanoengineering.These parameters cannot be measured experimentally because of nanocomposite type of materials [1].
There are many nondestructive material characterization techniques for obtaining permittivities of dielectric materials (see a short review of recent results in [2]).For the first time, the problem of complex permittivity determination in closed rectangular waveguides using propagation constant measurements was investigated in [3].For open dielectric waveguides of arbitrary cross-section such problems can be set up as inverse eigenvalue problems [4] of the theory of optical waveguides [5].Problems described in [4] are based on integral formulations of the underlying equations.
Inverse problems for determining of dielectric permittivity were studied in many works; see, for example, [6][7][8][9] and references therein.In [10], authors present reconstruction of complex permittivity using finite difference time-domain (FDTD) method.In [11,12], an analytical method for reconstruction of permittivity of a lossy arbitrary shaped body inside a three-dimensional waveguide using transmitted and reflected data was developed.Results of [11,12] are based on the volume singular integral equation (VSIE) method [8,9] for the system of Maxwell's equations.
The new approximately globally convergent method for reconstruction of permittivity function was developed in [13].This method was further verified on computationally simulated and on experimental data in [14][15][16][17][18][19][20][21][22].An adaptive version of the globally convergent method is developed and computationally verified for the first time in [23].In recent works [16,17,21], reconstruction of the spatially distributed dielectric functions and shapes of objects placed in the air as well as of buried objects in the dry sand from blind backscattered experimental data using two-stage numerical procedure of [13,15,24] is presented.In the first stage in [16,21,22], the approximately globally convergent method of [13] is applied to get a good first approximation for the exact permittivity function.Then, the local adaptive finite element method of [25,26] is used in the second stage to refine the solution obtained in the first stage.Results of this stage are presented in [15,17].An approximately globally convergent algorithm in a frequency domain for the case of reconstruction of dielectric function in an optical fiber is recently presented in [27].
The method investigated in this work significantly differs from the above described techniques.Compared with [16][17][18][19][20][21]27], we use measurements of the propagation constant 2 Physics Research International instead of the electric field data which is typical for waveguides applications.Our new algorithms use solution of integral equation over the boundary of a waveguide and thus they are less computationally expensive compared with methods of [11,12] which are also based on the volume singular integral equation (VSIE) technique [8,9].The models of [11,12] require the solution of three-dimensional vector problems and thus they are more computationally demanding and need powerful computer resources.
We note that two mathematical models of optical waveguides were investigated in detail by methods of integral equations in [19,20,[28][29][30][31][32][33][34].In [19,20,[28][29][30][31], the model of step-index waveguides was studied and, in [32][33][34], the model of waveguides without a sharp boundary was studied.A review of modern results in this field is given in [35].To solve our eigenvalue problem, we are not requiring the information about specific values of eigenfunctions.In our algorithms, it is enough only to know that the fundamental mode is excited and then to measure its propagation constant for one or two frequencies.Main applications of our algorithm are, for example, in detection of defects in metamaterials and in nanoelectronics [36,37], calculation of dielectric constant of saline water [38,39], and in microwave imaging technology in remote sensing [40].
An outline of this paper is a follows.In Section 2, we first present a mathematical model of a weakly guiding stepindex arbitrarily shaped optical fiber and formulate three inverse spectral problems.Then, we present new numerical algorithms for calculation of dielectric constants based on an approximate solution of a nonlinear nonselfadjoint eigenvalue problem for a system of weakly singular integral equations.In Section 3, we present reconstruction results which show efficiency and robustness of new developed algorithms.

Nonlinear Eigenvalue Problem for Transverse Wavenumbers.
Let us consider an optical fiber as a regular cylindrical dielectric waveguide in a free space.The cross-section of the waveguide's core is a bounded domain Ω  with twice continuously differentiable boundary  (see Figure 1).The axis of the cylinder is parallel to the  3 -axis.Let Ω  =  2 \ Ω  be the unbounded domain of the cladding.Let the permittivity be prescribed as a positive piecewise constant function  which is equal to a constant  ∞ in the domain Ω  and to a constant  + >  ∞ in the domain Ω  .
Eigenvalue problems of optical waveguide theory [5] are formulated on the base of the set of homogeneous Maxwell equations: Here, E and H are electric and magnetic field vectors and  0 and  0 are the free-space dielectric and magnetic constants.Nontrivial solutions of set (1), which have the form are called the eigenmodes of the waveguide.Here, positive  is the radian frequency,  is the propagation constant, E and H are complex amplitudes of E and H, and  = ( 1 ,  2 ).
In forward eigenvalue problems, the permittivity is known and it is necessary to calculate longitudinal wavenumbers  =  √  0  0 and propagation constants  such that there exist eigenmodes.The eigenmodes have to satisfy to a transparency condition at the boundary  and to a condition at infinity.
In inverse eigenvalue problems, it is necessary to reconstruct the unknown permittivity  by some information on natural eigenmodes which exist for some eigenvalues  and .The main question is how many observations of natural eigenmodes are enough for unique and stable reconstruction of the permittivity.
The domain Ω  is unbounded.Therefore, it is necessary to formulate a condition at infinity for complex amplitudes E and H of eigenmodes.Let us confine ourselves to the investigation of the surface modes only.The propagation constants  of surface modes are real and belong to the interval  = ( ∞ ,  + ).The amplitudes of surface modes satisfy the following condition: Here,  = √ 2 −  2  ∞ > 0 is the transverse wavenumber in the cladding.
Denote by  = √ 2  + −  2 the transverse wavenumber in the waveguide's core.Under the weakly guidance approximation [5], the original problem was reduced in [29] to the calculation of numbers  and  such that there exist nontrivial solutions  =  1 =  2 of Helmholtz equations: which satisfy the transparency conditions Here, ] is the normal derivative on ,  − (resp.,  + ) is the limit value of a function  from the interior (resp., the exterior) of .
Denote by  the space of continuous and continuously differentiable functions in Ω  and Ω  and twice continuously differentiable functions in Ω  and Ω  , satisfying to condition (3).Let us calculate nontrivial solutions  of problem ( 4)-( 5) in the space .
Let  > 0 be a given number.A nonzero function  ∈  is called an eigenfunction of problem ( 4)-( 5) corresponding to a real eigenvalue  if relationships (4)-( 5) hold.
The next theorem follows from results of [29].
Theorem 1.For any  > 0, the eigenvalues  of problem (4)-( 5) can be only positive isolated numbers.Each number  depends continuously on .
For waveguides of circular cross-section, the analogous results about the localization of the surface modes spectrum and about the continuous dependence between the transverse wavenumbers  and  were obtained in [5].The results of [5] are valid only for waveguides of circular cross-section by the method of separation of variables.Theorem 1 generalizes the results of [5] for the case of an arbitrary smooth boundary.The next theorem follows from results of [34] (see an illustration in Figure 2).Theorem 2. The following statements hold.
(2) For any  > 0, the smallest eigenvalue  1 () is positive and simple (its multiplicity is equal to one); corresponding eigenfunction  1 can be chosen as the positive function in the domain Ω  .
For a given  > 0, the smallest eigenvalue  1 () and corresponding eigenfunction  1 define the eigenmode which is called the fundamental mode (see the bottom curve plotted by the red solid line in Figure 2).Thus, Theorem 2 states, particularly, that, for any  > 0, there exists exactly one fundamental mode.
Let  > 0 be a given number.A nonzero element  ∈  is called an eigenfunction of the operator-valued function () corresponding to a characteristic value  > 0, if (15) holds.
It follows from results of [29] that original problem (4)-( 5) is spectrally equivalent to (15).Namely, the following theorem is true.Theorem 3. Let  > 0 be a given number.The following statements hold.
(1) If a function  is the eigenfunction of the operatorvalued function () corresponding to a characteristic value  0 > 0, then the function  defined by equalities (6), where  =  0 , is the eigenfunction of problem (4)-( 5) corresponding to the eigenvalue  0 .
(2) Each eigenfunction of problem (4)-( 5) corresponding to an eigenvalue  0 > 0 can be represented in the form of singlelayer potentials (6) with some Hölder-continuous densities  and , respectively.At the same time, the function
Let us formulate the nonlinear spectral problem for transverse wavenumbers as follows.Suppose that the boundary  of the waveguide's cross-section and the number  > 0 are given.It is necessary to calculate all characteristic values  of the operator-valued function () in some given interval.
A spline-collocation method was proposed in [31] for numerical calculations of characteristic values  of the operator-valued function () for given , and problem (15) for each fixed  was reduced to a nonlinear finite-dimensional eigenvalue problem of the form where  is the number of collocation points.For numerical solution of obtained nonlinear finite-dimensional eigenvalue problem, we use the residual inverse iteration method [42].Any iterative numerical method for computations of the nonlinear eigenvalues  needs good initial approximations for each given .Usually initial approximations are chosen by a physical intuition using prior information on solutions.If we model fundamentally new types of waveguides, solve inverse problems, or investigate defects in fibers, then we may not have accurate prior information on solutions.
In this case, we can investigate spectral properties of the matrix   () as a function of variable  for each fixed .For each given point in an investigated domain of parameters  and , we calculate the condition number of matrix   : where  1 and   are maximal and minimal singular values of matrix.If for given  a number  is equal to a nonlinear eigenvalue of   (), then the condition number is equal to infinity.Therefore, the numbers  for which the condition number is big enough are good approximations for eigenvalues.Singular values are calculated by singular value decomposition (SVD) method: where ,  are unitary matrices and  is a diagonal matrix; the singular numbers form . The calculations are based on unitary transformations of the matrix  and therefore are stable.
In the next step, for each  in the investigated interval, we use obtained initial approximations for numerical calculations of nonlinear eigenvalues  by the residual inverse iteration method.
On the base of solution of nonlinear eigenvalue problem (15) for transverse wavenumbers, we solve the forward and inverse spectral problems for weakly guiding step-index waveguides.

Forward Spectral Problem.
Clearly, if for given  the characteristic values  of the operator-valued function () were calculated and also if the permittivities  + ,  ∞ are known, then the longitudinal wavenumber  and the propagation constant  are calculated by the following explicit formulas: For each given  + ,  ∞ , and , the function   () generates a function  2 of variable  2 .As an example in Figure 3, we present the plot of the computed function  2 =  2 ( 2 ) corresponding to the fundamental mode (see the bottom curve plotted by the red solid line in Figure 2) of the circular waveguide.Here,  + = 2,  ∞ = 1, and  = 1.We will use two points marked by circle and by square in the next sections as test points for numerical solution of inverse spectral problems.

Inverse Spectral Problems.
In this subsection, we present three algorithms for approximate solution of three inverse spectral problems.The algorithms are based on previous numerical solution of nonlinear eigenvalue problem (15) for transverse wavenumbers, namely, on numerical calculations of characteristic values  1 () of the operator-valued function () for fundamental mode and for  in an appropriate interval.

Algorithm for Reconstruction of the Permittivity of the
Core.The first inverse problem is formulated as follows.Suppose that the boundary  of the waveguide's cross-section and the permittivity  ∞ of the cladding are given.Suppose that the propagation constant  of the fundamental mode is measured for a given wavenumber .The measuring can be done by experimental methods described, for example, in [3].It is necessary to calculate the permittivity  + of the waveguide's core.
The solution of this inverse spectral problem is calculated in the following way.First, we compute the number which is calculated for given values of , , and  ∞ .Then, the transverse wavenumber () is calculated by the splinecollocation method for the obtained  of the fundamental mode.This number is calculated by an interpolation of the function  1 () with respect to the points obtained when the nonlinear spectral problem for transverse wavenumbers was numerically solved.Finally, the permittivity of the waveguide's core is calculated by the following explicit formula:

Algorithm for Reconstruction of the Permittivity of the
Cladding.Suppose that the permittivity of the core is given and that the propagation constant  of the fundamental mode is measured for a given wavenumber .It is necessary to calculate the permittivity  ∞ of the waveguide's cladding.The unique exact solution of the inverse spectral problem is marked by the red circle.The approximate solution obtained by the splinecollocation method coincides with the red circle for the exact .
For randomly noised β, approximate solutions are intersections of dashed lines.They belong to the red rhomb.
The solution of this problem is calculated in the following way.First, we compute the number which is calculated for given values of , , and  + .Second, the transverse wavenumber  is calculated by the splinecollocation method for the obtained  for the fundamental mode.Third, the permittivity of the waveguide's cladding is calculated by the following explicit formula:

Algorithm for Reconstruction of Two Permittivities of the
Core and of the Cladding.The full variant of our problem is the reconstruction of both permittivities of the core and of the cladding.The solution for the fundamental mode of nonlinear eigenvalue problem (15) for transverse wavenumbers gives us an implicit function  + of the variable  ∞ for each fixed longitudinal wavenumber  and propagation constant .For example, in Figure 4, the blue and the green solid lines are plots of this function for given pairs of  and .The intersection of these lines marked by the red circle is the unique exact solution of our problem.Therefore, we calculate the permittivities  + and  ∞ as the solution of the following nonlinear system of two equations: Here,  1 and  2 are some given distinct longitudinal wavenumbers,  1 and  2 are corresponding measured propagation constants, and  is the function of variable  ∞ for fixed   and   ,  = 1, 2.

Results and Discussion
Let us describe numerical results obtained for a nonlinear eigenvalue problem (15) for transverse wavenumbers.We compare our numerical results with the well-known exact solution for the circular waveguide obtained by the method of separation of variables (see, e.g., [5]).In Figure 2, we present some dispersion curves for surface eigenwaves of the circular waveguide.Here, the exact solution is plotted by solid lines.We started our numerical calculations with computations of initial approximations for nonlinear eigenvalues .To do that, we have used SVD as was described in Section 2.1.Calculated initial approximations are marked in Figure 2 by blue circles.The blue circles are only the initial approximations for eigenvalues .We apply these initial approximations as start points in the residual inverse iteration method.Using this iteration method, for each given , we solved numerically the finite-dimensional nonlinear eigenvalue problem (20) depending on eigenvalues .The numerical solution which resulted in this method is marked in Figure 2 by red crosses.
Note here that, applying previously calculated initial approximations for nonlinear eigenvalues , we numerically solve this problem directly and without any prior physical information.Figure 2 illustrates that the initial approximations are good enough and we can use them in the majority of cases without iterative refinement of solutions.
Our next numerical experiment shows reconstruction of core's permittivity using the algorithm in Section 2.3.1.The mathematical analysis of existence of the solution of the original spectral problem (4)-( 5) for transverse wavenumbers is presented in Theorems 1 and 2.An illustration of the theoretical results is shown in Figure 2. Analyzing Figure 2, we observe that the fundamental mode (see the red solid curve in Figure 2) exists for each wavenumber  > 0. The fundamental mode is unique; its dispersion curve does not intersect with any other curves and well separated from them.Therefore, the inverse spectral problem's solution exists and is unique for each wavenumber  > 0; this solution depends continuously on given data.In other words, the inverse spectral problem is well posed by Hadamard.We can see an example of this continuity dependence in Figure 5.The red solid line is the plot of the function  + of  2 for the fundamental mode.In our computations, by analogy with [13,14], we have introduced random noise in the propagation constant as where  = 3.98518 is the exact measured propagation constant,  ∈ (−1, 1) are randomly distributed numbers, and  is the noise level.In our computations, we have used  = 0.05 and, thus, the noise level was 5%.Some numerical results of reconstruction of  + are presented in Figure 5.The approximated value of  + for the noise-free data is marked in Figure 5 by the blue circle.Approximated values of  + for randomly distributed noise β with the 5% noise level are marked in Figure 5 by the black points.Using this figure, we observe that the approximate solutions even for the randomly noised β were stable.We can conclude that, for the unique and stable reconstruction of the constant waveguide permittivity  + , it is enough to measure the propagation constant  of the fundamental mode for only one wavenumber .
Absolutely analogous results were obtained for reconstruction of cladding's permittivity using the algorithm in Section 2.3.2.As in Figure 5, the red solid line in Figure 6 is the plot of continuous function  ∞ of squared  for the fundamental mode.The approximate solution obtained by the spline-collocation method is marked by the blue circle for the exact  and by black points for the randomly noised β.Using this figure, we observe that the approximate solutions for the randomly noised β were stable in this case too.For  this test, we can also conclude that, for the unique and stable reconstruction of the constant permittivity  ∞ , it is enough to measure the propagation constant  of the fundamental mode for only one wavenumber .
Finally, we show simultaneous reconstruction of both permittivities of core and of cladding using the algorithm in Section 2.3.3.The approximate solution obtained by the spline-collocation method is marked in Figure 4 by the red circle for the exact propagation constant.We have introduced the random noise in the propagation constant similarly with above described tests.The approximate solutions for the noised β are presented in Figure 4 as intersections of dashed lines.Using Figure 4, we observe that the approximate solutions for the randomly noised  belong to the red rhomb and are stable.Therefore, we can conclude that, for the unique and stable reconstruction of the dielectric constants  ∞ and  + , it is enough to measure the propagation constant  of the fundamental mode for two distinct wavenumbers .

Conclusion
In this work, we have formulated three inverse spectral problems and proved that they are well posed.It is important to note that, in our algorithms, any information on specific values of eigenfunctions is not required.For solution of these inverse problems, it is enough to know that the fundamental mode is excited and then to measure its propagation constant for one or for two frequencies.This approach satisfies

Figure 1 :
Figure 1: The cross-section of a cylindrical dielectric waveguide in a free space.

Figure 2 :
Figure 2: The plot of the dispersion curves computed by the spline-collocation method for surface eigenmodes in a weakly guiding dielectric waveguide of the circular cross-section.The exact solutions are plotted by solid lines.The solution for the fundamental mode is plotted by the lower red solid line.Numerical solutions of the residual inverse iteration method are marked by red crosses.Initial approximations for SVDs are marked by blue circles.

Figure 3 :
Figure 3: The red solid line presents the plot of the computed function  2 =  2 ( 2 ) corresponding to the fundamental mode of the circular waveguide.

Figure 4 :
Figure 4: The blue and the green solid lines are plots of function  + =  + ( ∞ ) for given pairs of  and  of the fundamental mode.The unique exact solution of the inverse spectral problem is marked by the red circle.The approximate solution obtained by the splinecollocation method coincides with the red circle for the exact .For randomly noised β, approximate solutions are intersections of dashed lines.They belong to the red rhomb.

Figure 5 :
Figure 5: The red solid line is the plot of function  + =  + ( 2 ) for the fundamental mode.The approximate solution obtained by the spline-collocation method is marked by the blue circle for the exact  and by black points for the randomly noised β.

Figure 6 :
Figure 6: The red solid line is the plot of function  ∞ =  ∞ ( 2 ) for the fundamental mode.The approximate solution obtained by the spline-collocation method is marked by the blue circle for the exact  and by black points for the randomly noised β.