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

Authors: E. M. Karchevskii; A. O. Spiridonov; A. I. Repina; L. Beilina
Journal: Physics Research International
 (2014)

---

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


---

## Body

## 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 Section2, 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Ωi with twice continuously differentiable boundary γ (see Figure 1). The axis of the cylinder is parallel to the x3-axis. Let Ωe=R2∖Ω¯i 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 Ωe and to a constant ε+>ε∞ in the domain Ωi.Figure 1
The cross-section of a cylindrical dielectric waveguide in a free space.Eigenvalue problems of optical waveguide theory [5] are formulated on the base of the set of homogeneous Maxwell equations:
(1)rot⁡E=-μ0∂H∂t,rot⁡H=ε0ɛ∂E∂t.
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
(2)[EH](x,x3,t)=Re([EH](x)ei(βx3-ωt)),
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 x=(x1,x2).In forward eigenvalue problems, the permittivity is known and it is necessary to calculate longitudinal wavenumbersk=ωε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 k and β. The main question is how many observations of natural eigenmodes are enough for unique and stable reconstruction of the permittivity.The domainΩe 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 G=(kε∞,kε+). The amplitudes of surface modes satisfy the following condition:
(3)[EH]=e-σrO(1r),r=|x|⟶∞.
Here, σ=β2-k2ε∞>0 is the transverse wavenumber in the cladding.Denote byχ=k2ε+-β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 u=H1=H2 of Helmholtz equations:
(4)Δu+χ2u=0,x∈Ωi,Δu-σ2u=0,x∈Ωe,
which satisfy the transparency conditions
(5)u+=u-,∂u+∂ν=∂u-∂ν,x∈γ.
Here, ν is the normal derivative on γ, u- (resp., u+) is the limit value of a function u from the interior (resp., the exterior) of γ.Denote byU the space of continuous and continuously differentiable functions in Ωi and Ωe and twice continuously differentiable functions in Ωi and Ωe, satisfying to condition (3). Let us calculate nontrivial solutions u of problem (4)-(5) in the space U.Letσ>0 be a given number. A nonzero function u∈U 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).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.Theorem 2.
The following statements hold.
(1) For anyσ>0, there exist the denumerable set of positive eigenvalues χi(σ), where i=1,2,…, of a finite multiplicity with only cumulative point at infinity.
(2) For anyσ>0, the smallest eigenvalue χ1(σ) is positive and simple (its multiplicity is equal to one); corresponding eigenfunction u1 can be chosen as the positive function in the domain Ωi.
(3)χ1(σ)→0 at σ→0.For a givenσ>0, the smallest eigenvalue χ1(σ) and corresponding eigenfunction u1 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.To compute eigenmodes, we use the representation of eigenfunctions of problem (4)-(5) in the form of single-layer potentials [41]:
(6)u(x)=∫γΦ(χ;x,y)f(y)dl(y),x∈Ωi,u(x)=∫γΨ(σ;x,y)g(y)dl(y),x∈Ωe.
Here,
(7)Φ(χ;x,y)=i4H0(1)(χ|x-y|),Ψ(σ;x,y)=12πK0(σ|x-y|)
are the fundamental solutions of Helmholtz equations (4) correspondingly, H0(1)(z) is Hankel function of the first kind, K0(z) is Macdonald function, and unknown densities f and g belong to the Hölder space C0,α(γ). Then, f and g satisfy to the following system of boundary integral equations for x∈γ:
(8)A11(χ)f+A12(σ)g=0,A21(χ)f+A22(σ)g=0,
where
(9)(A11(χ)f)(x)=∫γΦ(χ;x,y)f(y)dl(y),(A12(σ)g)(x)=-∫γΨ(σ;x,y)g(y)dl(y),(A21(χ)f)(x)=12f(x)+∫γ∂Φ(χ;x,y)∂ν(x)f(y)dl(y),(A22(σ)g)(x)=12g(x)-∫γ∂Ψ(σ;x,y)∂ν(x)g(y)dl(y).Suppose that the boundaryγ is parametrically defined by a function r=r(t) of t∈[0,2π] and this parametrization is regular. We consider functions from C0,α(γ) and C1,α(γ) as Hölder-continues and Hölder-continuously differentiable 2π-periodic functions of parameter t. Then, for x=x(t)∈γ, where t∈[0,2π], system (8) is transformed to the following system of equations:
(10)Lp1+B11(χ,σ)p1+B12(χ,σ)p2=0,p2+B21(χ,σ)p1+B22(χ,σ)p2=0.
Here, new unknown functions are defined as follows:
(11)p1(τ)=(f(y)-g(y))|r′(τ)|,p2(τ)=f(y)+g(y),y=y(τ),τ∈[0,2π].
Integral operators are defined by the following relationships for t∈[0,2π]:
(12)(Lp1)(t)=-12π∫02πln⁡|sint-τ2|p1(τ)dτ,(Bij(χ,σ)p(j))(t)=12π∫02πhij(χ,σ;t,τ)p(j)(τ)dτ,
where
(13)h11(χ,σ;t,τ)=2π(G11(χ;t,τ)+G12(σ;t,τ)),h12(χ,σ;t,τ)=2π(G11(χ;t,τ)-G12(σ;t,τ))|r′(τ)|,h21(χ,σ;t,τ)=4π(G21(χ;t,τ)+G22(σ;t,τ)),h22(χ,σ;t,τ)=4π(G21(χ;t,τ)-G22(σ;t,τ))|r′(τ)|,G11(χ;t,τ)=Φ(χ;x,y)+12πln⁡|sint-τ2|,G12(σ;t,τ)=Ψ(σ;x,y)+12πln⁡|sint-τ2|,G21(χ;t,τ)=∂Φ(χ;x,y)∂ν(x),G22(σ;t,τ)=∂Ψ(σ;x,y)∂ν(x).
The linear operator L:C0,α(γ)→C1,α(γ) is continuous and continuously invertible (for details, see, e.g., [29]). It was proved in [29] that, for each χ>0,  σ>0, the linear operators
(14)B21(χ,σ),B22(χ,σ):C0,α(γ)⟶C0,α(γ),B11(χ,σ),B12(χ,σ):C0,α(γ)⟶C1,α(γ)
are compact. Therefore, system (8) is equivalent to the following Fredholm operator equation of the second kind:
(15)A(χ,σ)w=(I+B(χ,σ))w=0,
where w=(w1,w2)T,
(16)w1=Lp1∈C1,α(γ),w2=p2∈C0,α(γ),
and the compact operator B acting in the Banach space W=C1,α(γ)×C0,α(γ) is defined as follows:
(17)Bw=[B11L-1B12B21L-1B22][w1w2];
by I the identical operator is denoted.Letσ>0 be a given number. A nonzero element w∈W is called an eigenfunction of the operator-valued function A(χ) 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 functionw is the eigenfunction of the operator-valued function A(χ) corresponding to a characteristic value χ0>0, then the function u defined by equalities (6), where χ=χ0,
(18)f=w22+L-1w12|r′|,g=w22-L-1w12|r′|,
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 single-layer potentials (6) with some Hölder-continuous densities f and g, respectively. At the same time, the function
(19)w=(L((f-g)|r′|),f+g)
is the eigenfunction of the operator-valued function A(χ) corresponding to the characteristic value χ0.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 A(χ) in some given interval.A spline-collocation method was proposed in [31] for numerical calculations of characteristic values χ of the operator-valued function A(χ) for given σ, and problem (15) for each fixed σ was reduced to a nonlinear finite-dimensional eigenvalue problem of the form
(20)An(χ)wn=0,
where n 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 matrixAn(χ) 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 An:
(21)cond(An(χ))=ρ1ρn,
where ρ1 and ρn are maximal and minimal singular values of matrix. If for given σ a number χ is equal to a nonlinear eigenvalue of An(χ), 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:
(22)An(χ)=USV,S=diag⁡(ρ1,ρ2,…,ρn),
where U, V are unitary matrices and S is a diagonal matrix; the singular numbers form S. The calculations are based on unitary transformations of the matrix A 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 A(χ) were calculated and also if the permittivities ε+, ε∞ are known, then the longitudinal wavenumber k and the propagation constant β are calculated by the following explicit formulas:
(23)k2=σ2+χ2ε+-ε∞,β2=ε+σ2+ε∞χ2ε+-ε∞.
For each given ε+, ε∞, and i, the function χi(σ) generates a function β2 of variable k2. As an example in Figure 3, we present the plot of the computed function β2=β2(k2) 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 i=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.Figure 3
The red solid line presents the plot of the computed functionβ2=β2(k2) corresponding to the fundamental mode of the circular waveguide.


### 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 χ1(σ) of the operator-valued function A(χ) 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 k. 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(24)σ=β2-k2ε∞,
which is calculated for given values of β, k, 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 χ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:
(25)ε+=χ2+β2k2.


#### 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 k. 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(26)χ=k2ε+-β2,
which is calculated for given values of β, k, 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:
(27)ε∞=β2-σ2k2.


#### 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 k and propagation constant β. For example, in Figure 4, the blue and the green solid lines are plots of this function for given pairs of k and β.Figure 4
The blue and the green solid lines are plots of functionε+=ε+(ε∞) for given pairs of k 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 spline-collocation 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 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:
(28)χ(β12-k12ε∞)=k12ε+-β12,χ(β22-k22ε∞)=k22ε+-β22.
Here, k1 and k2 are some given distinct longitudinal wavenumbers, β1 and β2 are corresponding measured propagation constants, and χ is the function of variable ε∞ for fixed kj and βj,  j=1,2.


## 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Ωi with twice continuously differentiable boundary γ (see Figure 1). The axis of the cylinder is parallel to the x3-axis. Let Ωe=R2∖Ω¯i 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 Ωe and to a constant ε+>ε∞ in the domain Ωi.Figure 1
The cross-section of a cylindrical dielectric waveguide in a free space.Eigenvalue problems of optical waveguide theory [5] are formulated on the base of the set of homogeneous Maxwell equations:
(1)rot⁡E=-μ0∂H∂t,rot⁡H=ε0ɛ∂E∂t.
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
(2)[EH](x,x3,t)=Re([EH](x)ei(βx3-ωt)),
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 x=(x1,x2).In forward eigenvalue problems, the permittivity is known and it is necessary to calculate longitudinal wavenumbersk=ωε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 k and β. The main question is how many observations of natural eigenmodes are enough for unique and stable reconstruction of the permittivity.The domainΩe 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 G=(kε∞,kε+). The amplitudes of surface modes satisfy the following condition:
(3)[EH]=e-σrO(1r),r=|x|⟶∞.
Here, σ=β2-k2ε∞>0 is the transverse wavenumber in the cladding.Denote byχ=k2ε+-β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 u=H1=H2 of Helmholtz equations:
(4)Δu+χ2u=0,x∈Ωi,Δu-σ2u=0,x∈Ωe,
which satisfy the transparency conditions
(5)u+=u-,∂u+∂ν=∂u-∂ν,x∈γ.
Here, ν is the normal derivative on γ, u- (resp., u+) is the limit value of a function u from the interior (resp., the exterior) of γ.Denote byU the space of continuous and continuously differentiable functions in Ωi and Ωe and twice continuously differentiable functions in Ωi and Ωe, satisfying to condition (3). Let us calculate nontrivial solutions u of problem (4)-(5) in the space U.Letσ>0 be a given number. A nonzero function u∈U 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).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.Theorem 2.
The following statements hold.
(1) For anyσ>0, there exist the denumerable set of positive eigenvalues χi(σ), where i=1,2,…, of a finite multiplicity with only cumulative point at infinity.
(2) For anyσ>0, the smallest eigenvalue χ1(σ) is positive and simple (its multiplicity is equal to one); corresponding eigenfunction u1 can be chosen as the positive function in the domain Ωi.
(3)χ1(σ)→0 at σ→0.For a givenσ>0, the smallest eigenvalue χ1(σ) and corresponding eigenfunction u1 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.To compute eigenmodes, we use the representation of eigenfunctions of problem (4)-(5) in the form of single-layer potentials [41]:
(6)u(x)=∫γΦ(χ;x,y)f(y)dl(y),x∈Ωi,u(x)=∫γΨ(σ;x,y)g(y)dl(y),x∈Ωe.
Here,
(7)Φ(χ;x,y)=i4H0(1)(χ|x-y|),Ψ(σ;x,y)=12πK0(σ|x-y|)
are the fundamental solutions of Helmholtz equations (4) correspondingly, H0(1)(z) is Hankel function of the first kind, K0(z) is Macdonald function, and unknown densities f and g belong to the Hölder space C0,α(γ). Then, f and g satisfy to the following system of boundary integral equations for x∈γ:
(8)A11(χ)f+A12(σ)g=0,A21(χ)f+A22(σ)g=0,
where
(9)(A11(χ)f)(x)=∫γΦ(χ;x,y)f(y)dl(y),(A12(σ)g)(x)=-∫γΨ(σ;x,y)g(y)dl(y),(A21(χ)f)(x)=12f(x)+∫γ∂Φ(χ;x,y)∂ν(x)f(y)dl(y),(A22(σ)g)(x)=12g(x)-∫γ∂Ψ(σ;x,y)∂ν(x)g(y)dl(y).Suppose that the boundaryγ is parametrically defined by a function r=r(t) of t∈[0,2π] and this parametrization is regular. We consider functions from C0,α(γ) and C1,α(γ) as Hölder-continues and Hölder-continuously differentiable 2π-periodic functions of parameter t. Then, for x=x(t)∈γ, where t∈[0,2π], system (8) is transformed to the following system of equations:
(10)Lp1+B11(χ,σ)p1+B12(χ,σ)p2=0,p2+B21(χ,σ)p1+B22(χ,σ)p2=0.
Here, new unknown functions are defined as follows:
(11)p1(τ)=(f(y)-g(y))|r′(τ)|,p2(τ)=f(y)+g(y),y=y(τ),τ∈[0,2π].
Integral operators are defined by the following relationships for t∈[0,2π]:
(12)(Lp1)(t)=-12π∫02πln⁡|sint-τ2|p1(τ)dτ,(Bij(χ,σ)p(j))(t)=12π∫02πhij(χ,σ;t,τ)p(j)(τ)dτ,
where
(13)h11(χ,σ;t,τ)=2π(G11(χ;t,τ)+G12(σ;t,τ)),h12(χ,σ;t,τ)=2π(G11(χ;t,τ)-G12(σ;t,τ))|r′(τ)|,h21(χ,σ;t,τ)=4π(G21(χ;t,τ)+G22(σ;t,τ)),h22(χ,σ;t,τ)=4π(G21(χ;t,τ)-G22(σ;t,τ))|r′(τ)|,G11(χ;t,τ)=Φ(χ;x,y)+12πln⁡|sint-τ2|,G12(σ;t,τ)=Ψ(σ;x,y)+12πln⁡|sint-τ2|,G21(χ;t,τ)=∂Φ(χ;x,y)∂ν(x),G22(σ;t,τ)=∂Ψ(σ;x,y)∂ν(x).
The linear operator L:C0,α(γ)→C1,α(γ) is continuous and continuously invertible (for details, see, e.g., [29]). It was proved in [29] that, for each χ>0,  σ>0, the linear operators
(14)B21(χ,σ),B22(χ,σ):C0,α(γ)⟶C0,α(γ),B11(χ,σ),B12(χ,σ):C0,α(γ)⟶C1,α(γ)
are compact. Therefore, system (8) is equivalent to the following Fredholm operator equation of the second kind:
(15)A(χ,σ)w=(I+B(χ,σ))w=0,
where w=(w1,w2)T,
(16)w1=Lp1∈C1,α(γ),w2=p2∈C0,α(γ),
and the compact operator B acting in the Banach space W=C1,α(γ)×C0,α(γ) is defined as follows:
(17)Bw=[B11L-1B12B21L-1B22][w1w2];
by I the identical operator is denoted.Letσ>0 be a given number. A nonzero element w∈W is called an eigenfunction of the operator-valued function A(χ) 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 functionw is the eigenfunction of the operator-valued function A(χ) corresponding to a characteristic value χ0>0, then the function u defined by equalities (6), where χ=χ0,
(18)f=w22+L-1w12|r′|,g=w22-L-1w12|r′|,
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 single-layer potentials (6) with some Hölder-continuous densities f and g, respectively. At the same time, the function
(19)w=(L((f-g)|r′|),f+g)
is the eigenfunction of the operator-valued function A(χ) corresponding to the characteristic value χ0.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 A(χ) in some given interval.A spline-collocation method was proposed in [31] for numerical calculations of characteristic values χ of the operator-valued function A(χ) for given σ, and problem (15) for each fixed σ was reduced to a nonlinear finite-dimensional eigenvalue problem of the form
(20)An(χ)wn=0,
where n 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 matrixAn(χ) 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 An:
(21)cond(An(χ))=ρ1ρn,
where ρ1 and ρn are maximal and minimal singular values of matrix. If for given σ a number χ is equal to a nonlinear eigenvalue of An(χ), 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:
(22)An(χ)=USV,S=diag⁡(ρ1,ρ2,…,ρn),
where U, V are unitary matrices and S is a diagonal matrix; the singular numbers form S. The calculations are based on unitary transformations of the matrix A 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 A(χ) were calculated and also if the permittivities ε+, ε∞ are known, then the longitudinal wavenumber k and the propagation constant β are calculated by the following explicit formulas:
(23)k2=σ2+χ2ε+-ε∞,β2=ε+σ2+ε∞χ2ε+-ε∞.
For each given ε+, ε∞, and i, the function χi(σ) generates a function β2 of variable k2. As an example in Figure 3, we present the plot of the computed function β2=β2(k2) 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 i=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.Figure 3
The red solid line presents the plot of the computed functionβ2=β2(k2) corresponding to the fundamental mode of the circular waveguide.


## 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 χ1(σ) of the operator-valued function A(χ) 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 k. 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(24)σ=β2-k2ε∞,
which is calculated for given values of β, k, 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 χ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:
(25)ε+=χ2+β2k2.


### 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 k. 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(26)χ=k2ε+-β2,
which is calculated for given values of β, k, 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:
(27)ε∞=β2-σ2k2.


### 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 k and propagation constant β. For example, in Figure 4, the blue and the green solid lines are plots of this function for given pairs of k and β.Figure 4
The blue and the green solid lines are plots of functionε+=ε+(ε∞) for given pairs of k 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 spline-collocation 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 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:
(28)χ(β12-k12ε∞)=k12ε+-β12,χ(β22-k22ε∞)=k22ε+-β22.
Here, k1 and k2 are some given distinct longitudinal wavenumbers, β1 and β2 are corresponding measured propagation constants, and χ is the function of variable ε∞ for fixed kj and βj,  j=1,2.


## 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 k. 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(24)σ=β2-k2ε∞,
which is calculated for given values of β, k, 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 χ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:
(25)ε+=χ2+β2k2.


## 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 k. 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(26)χ=k2ε+-β2,
which is calculated for given values of β, k, 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:
(27)ε∞=β2-σ2k2.


## 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 k and propagation constant β. For example, in Figure 4, the blue and the green solid lines are plots of this function for given pairs of k and β.Figure 4
The blue and the green solid lines are plots of functionε+=ε+(ε∞) for given pairs of k 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 spline-collocation 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 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:
(28)χ(β12-k12ε∞)=k12ε+-β12,χ(β22-k22ε∞)=k22ε+-β22.
Here, k1 and k2 are some given distinct longitudinal wavenumbers, β1 and β2 are corresponding measured propagation constants, and χ is the function of variable ε∞ for fixed kj and βj,  j=1,2.


## 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 Section2.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 k>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 k>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.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 β~.In our computations, by analogy with [13, 14], we have introduced random noise in the propagation constant as
(29)β~=β(1+pα),
where β=3.98518 is the exact measured propagation constant, α∈(-1,1) are randomly distributed numbers, and p is the noise level. In our computations, we have used p=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 k.Absolutely analogous results were obtained for reconstruction of cladding’s permittivity using the algorithm in Section2.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 k.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 β~.Finally, we show simultaneous reconstruction of both permittivities of core and of cladding using the algorithm in Section2.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 k.


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



---

*Source: 253435-2014-09-18.xml*