#### Abstract

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.

#### 1. 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–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–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–21, 27], we use measurements of the propagation constant 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–34]. In [19, 20, 28–31], the model of step-index waveguides was studied and, in [32–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 step-index 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.

#### 2. Methods of Solution

##### 2.1. 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 -axis. Let 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, and are electric and magnetic field vectors and and 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, and are complex amplitudes of and , and .

In forward eigenvalue problems, the permittivity is known and it is necessary to calculate longitudinal wavenumbers 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 and 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, is the transverse wavenumber in the cladding.

Denote by 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 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 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 , 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.**(1) For any , there exist the denumerable set of positive eigenvalues , where , of a finite multiplicity with only cumulative point at infinity.**(2) For any , the smallest eigenvalue is positive and simple (its multiplicity is equal to one); corresponding eigenfunction can be chosen as the positive function in the domain .**(3) at .*

For a given , the smallest eigenvalue and corresponding eigenfunction 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 , there exists exactly one fundamental mode.

To compute eigenmodes, we use the representation of eigenfunctions of problem (4)-(5) in the form of single-layer potentials [41]: Here, are the fundamental solutions of Helmholtz equations (4) correspondingly, is Hankel function of the first kind, is Macdonald function, and unknown densities and belong to the Hölder space . Then, and satisfy to the following system of boundary integral equations for : where

Suppose that the boundary is parametrically defined by a function of and this parametrization is regular. We consider functions from and as Hölder-continues and Hölder-continuously differentiable -periodic functions of parameter . Then, for , where , system (8) is transformed to the following system of equations: Here, new unknown functions are defined as follows: Integral operators are defined by the following relationships for : where The linear operator is continuous and continuously invertible (for details, see, e.g., [29]). It was proved in [29] that, for each , , the linear operators are compact. Therefore, system (8) is equivalent to the following Fredholm operator equation of the second kind: where , and the compact operator acting in the Banach space is defined as follows: by the identical operator is denoted.

Let be a given number. A nonzero element is called an eigenfunction of the operator-valued function corresponding to a characteristic value , 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 be a given number. The following statements hold.**(1) If a function is the eigenfunction of the operator-valued function corresponding to a characteristic value , then the function defined by equalities (6), where ,
**
is the eigenfunction of problem (4)-(5) corresponding to the eigenvalue .**(2) Each eigenfunction of problem (4)-(5) corresponding to an eigenvalue can be represented in the form of single-layer potentials (6) with some Hölder-continuous densities and , respectively. At the same time, the function
**
is the eigenfunction of the operator-valued function corresponding to the characteristic value .*

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 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 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.

##### 2.2. 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 of variable . As an example in Figure 3, we present the plot of the computed function corresponding to the fundamental mode (see the bottom curve plotted by the red solid line in Figure 2) of the circular waveguide. Here, , , and . 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.

##### 2.3. 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 of the operator-valued function for fundamental mode and for in an appropriate interval.

###### 2.3.1. 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 spline-collocation method for the obtained of the fundamental mode. This number is calculated by an interpolation of the function 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:

###### 2.3.2. 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 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 spline-collocation method for the obtained for the fundamental mode. Third, the permittivity of the waveguide’s cladding is calculated by the following explicit formula:

###### 2.3.3. 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, and are some given distinct longitudinal wavenumbers, and are corresponding measured propagation constants, and is the function of variable for fixed and , .

#### 3. 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 . 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 ; 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 for the fundamental mode.

In our computations, by analogy with [13, 14], we have introduced random noise in the propagation constant as where is the exact measured propagation constant, are randomly distributed numbers, and is the noise level. In our computations, we have used 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 .

#### 4. 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 the practice of physical experiments because usually the fundamental mode is excited in waveguides in real applications [3]. Moreover, the fundamental mode can be excited only for enough wide range of frequencies.

For the approximate solution of the inverse problems, we propose to solve the nonlinear spectral problem for transverse wavenumbers in order to compute the dispersion curve for the fundamental mode. These calculations are done accurately by the spline-collocation method. Finally, we have demonstrated unique and stable reconstruction of the permittivities using new inverse algorithms.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The authors would like to thank Michael Havrilla and George Hanson for helpful discussions concerning the topics covered in this work. This work was funded by the subsidy allocated to Kazan Federal University for the state assignment in the sphere of scientific activities; the work was supported also by RFBR and by Government of Republic Tatarstan, Grant 12-01-97012-r_povolzh’e_a. The research of Larisa Beilina was supported by the Swedish Research Council, the Swedish Foundation for Strategic Research (SSF) through the Gothenburg Mathematical Modelling Centre (GMMC), and by the Swedish Institute.