# Nonlinear Double-Layer Bragg Waveguide: Analytical and Numerical Approaches to Investigate Waveguiding Problem

Authors: Yury G. Smirnov; Eugenii Yu. Smol’kin; Dmitry V. Valovik
Journal: Advances in Numerical Analysis
 (2014)

---

## Abstract

The paper is concerned with propagation of surface TE waves in a circular nonhomogeneous two-layered dielectric waveguide filled with nonlinear medium. The problem is reduced to the analysis of a
nonlinear integral equation with a kernel in the form of the Green function. The existence of propagating TE waves for chosen nonlinearity (the Kerr law) is proved using the contraction mapping method. Conditions under whichk waves can propagate are obtained, and intervals of localization of the corresponding propagation constants are found. For numerical solution of the problem, a method based on solving an auxiliary Cauchy problem (the shooting method) is proposed. In numerical experiment, two types of nonlinearities are considered and compared: the Kerr nonlinearity and nonlinearity with saturation. New propagation regime is found.


---

## Body

## 1. Introduction

Theory of circle cylindrical dielectric waveguides attracts attention for a long time. Linear theory of such waveguides is known for years; see, for example, [1–4]. At the same time it is well known that the permittivity of a dielectric, in general, depends nonlinearly on the intensity of an electromagnetic field; see [5, 6]. For this reason, the linear theory can be applied only for fields of low intensity. However, for applications, sometimes it is necessary to raise the intensity of the field, for example, in order to compensate the losses. What happens in the case of high intensity (when the permittivity of the dielectric depends nonlinearly on the intensity of the field)? Is it possible to preserve waveguide regimes and, if so, how to determine propagation constants in the nonlinear case? It is not always easy to answer these questions. However, in the case of “simple” geometry (circle cylindrical and plane-layered waveguides) and polarized (TE and TM) electromagnetic waves, it is possible to answer these questions.To the best of our knowledge, the first rigorous study of polarized electromagnetic wave propagation in a nonlinear circle cylindrical dielectric homogeneous waveguide is in [7, 8]. Then, there were several works, where some important cases for nonlinear but homogeneous permittivity have been investigated; see [9–12]. The next step was to apply earlier developed technique to the cases of multilayered waveguides and inhomogeneous nonlinear permittivity.In the paper [13], we considered integral equation approach to derive dispersion equations in a nonlinear waveguiding problem. This approach can be used for numerical implementation. However, there exists a faster and simpler numerical approach that allows determining propagation constants in a wide range of waveguiding problems. Here, we are going to demonstrate this numerical approach for special waveguiding problem. As Bragg (or multilayered) waveguides have application (see, e.g., [14, 15]), we demonstrate here that the approach from [13] can be easily extended to be applied to more complicated problems. In order to justify the analytical approach given here, we will widely use the results of the paper [13].


## 2. Statement of the Problem

Let us consider a three-dimensional spaceℝ3 with cylindrical coordinate system Oρφz. The space is filled with an isotropic medium of dielectric permittivity ε3ε0=const, where ε0 is the permittivity of free space, without sources. The medium is assumed to be isotropic and nonmagnetic. A cylindrical dielectric waveguide of circular cross section is presented as
(1)W∶={(ρ,φ,z)∶0⩽ρ<R1,0⩽φ<2π}∪{(ρ,φ,z)∶R1⩽ρ⩽R2,0⩽φ<2π}and is placed inℝ3 with its generating line parallel to the axis Oz. It is supposed that everywhere μ=μ0, where μ0 is the permeability of free space.The cross section of the waveguide, which is perpendicular to its axis, consists of two concentric circles of radiiR1 and R2; that is, the waveguide is two-layered (see Figure 1): R1 is the radius of the inner cylinder and R2-R1 is the thickness of the outer cylindrical shell. In the case of linear homogeneous media, such multilayered waveguides were studied in [3]. A practical application of a nonlinear Bragg waveguide is pointed out in [15].Figure 1
Geometry of the problem.The geometry of the problem is shown in Figure1. The waveguide is unlimitedly continued in z direction.LetE and H be complex amplitudes of an electromagnetic field. The complex amplitudes E, H must satisfy Maxwell’s equations [7]
(2)rot⁡H=-iωɛE,rot⁡E=iωμH,
have continuous tangential field components on the media interface ρ=R1, ρ=R2 and obey the radiation condition at infinity; that is, the electromagnetic field decays exponentially as ρ→∞ in the region ρ>R2; ω is the circular frequency.The permittivity in the entire space has the formɛ=ɛ~ε0, where
(3)ɛ~={ε1,0⩽ρ<R1,ε2(ρ)+α~|E|2,R1⩽ρ⩽R2,ε3,ρ>R2,and|E|2=|(Ee-iωt,eρ)|2+|(Ee-iωt,eφ)|2+|(Ee-iωt,ez)|2; ek is the orthonormal vector in the k direction;  (·,·) is the Euclidean scalar product; ε2(ρ)∈C[R1,R2]; α~ is a real constant; ε1⩾1, ε3⩾1 are real positive constants such that
(4)max⁡(ε1,ε3)<ε~2=min⁡ρ∈[R1,R2]ε2(ρ).


## 3. TE Waves

Let us consider TE-polarized waves in the monochromatic mode (see [16]):
(5)Ee-iωt=e-iωt(0,Eφ,0)T,He-iωt=e-iωt(Hρ,0,Hz)T,
where E,H are complex amplitudes and
(6)Eφ=Eφ(ρ)eiγz,Hρ=Hρ(ρ)eiγz,Hz=Hz(ρ)eiγz,
where γ is unknown spectral parameter.Letk02∶=ω2με0. Thus, substituting components (6) into (2) and using the notation u(ρ;γ)∶=Eφ(ρ), we obtain
(7)(ρ-1(ρu)′)′+(k02ɛ~-γ2)u=0,
where ɛ~=ε0-1ɛ and ɛ~ is defined by formula (3).We assume that functionu is sufficiently smooth:
(8)u(ρ)∈C1[0,+∞)∩C2(0,R1)∩C2(R1,R2)∩C2(R2,+∞).We will seekγ under condition  γ2>max⁡(ε1,ε3).Letk12∶=γ2-k02ε1 and k32∶=γ2-k02ε3. In the domains ρ<R1 and ρ>R2, equation (7) takes the forms, respectively,
(9)u′′+ρ-1u′-ρ-2u-k12u=0,(10)u′′+ρ-1u′-ρ-2u-k32u=0.In the domainR1⩽ρ⩽R2, equation (7) takes the form
(11)u′′+ρ-1u′-ρ-2u+k2(ρ)u+αu3=0,
where α∶=α~ω2μ0ε0, k2(ρ)∶=k22(ρ)-γ2, k22(ρ)∶=k02ε2(ρ).The necessary solutions to (9) and (10) must be written in the following form:
(12)u=C1I1(k1ρ),ρ<R1,(13)u=C2K1(k3ρ),ρ>R2.The functionsI1 and K1 are the modified Bessel functions and C1 and C2 are constants. The radiation condition is fulfilled because K1(k3ρ)→0 exponentially as ρ→∞ (see [17]).


## 4. Transmission Conditions and Transmission Problem

Transmission conditions for the functionsu and u′ result from the continuity conditions for the tangential field components and have the form
(14)[u]|ρ=R1=0,[u′]|ρ=R1=0,[u]|ρ=R2=0,[u′]|ρ=R2=0,
where [v]|ρ=s=lim⁡ρ→s-0⁡(ρ)-lim⁡ρ→s+0⁡v(ρ) is the jump of the limit values of the function at a point s.Let us formulate the transmission eigenvalue problem (problemP) to which the problem of surface waves propagating in a cylindrical waveguide has been reduced. The goal is to find quantities γ such that, for given C1≠0 (or C2≠0), there is a nonzero function u(ρ;γ) that is defined by formulas (12) and (13) for ρ<R1 and ρ>R2, respectively, and solves equation (11) for R1<ρ<R2; moreover, the function u(ρ;γ) thus defined satisfies conditions (8) and transmission conditions (14).The quantities solving problemP are called eigenvalues, and the corresponding functions u(ρ;γ) are called eigenfunctions. Such definition of an eigenvalue was given in [18] and found application later in similar nonlinear problems (see, e.g., [19, 20]). It should be noted that the eigenvalue γ depends on the value of the eigenfunction on one of the waveguide boundaries.


## 5. Nonlinear Integral Equation

Let us consider the equation (linear part of (11))
(15)Lu=0,
where L=(d/dρ)(ρ(d/dρ))+(k2(ρ)ρ-ρ-1), ρ∈[R1,R2].Letλ∶=γ2. Consider the boundary value problem
(16)ρv′′+v′+(k22(ρ)ρ-ρ-1)v=λρv,ρ∈[R1,R2],v′|ρ=R1=v′|ρ=R2=0.Letλn,vn(ρ) be a complete system of (real) eigenvalues and orthonormalized with weight ρ eigenfunctions of boundary value problem (16) (see [21]); this system exists as coefficients of the equation satisfy all necessary conditions. Then, for λ≠λn, the boundary value problem
(17)Lu=0,u′|ρ=R1=u′|ρ=R2=0
has only the trivial solution. This means that, for λ≠λn, there uniquely exists the Green function G(ρ,s;λ) of the boundary value problem
(18)LG=-δ(ρ-s),∂ρG|ρ=R1=∂ρG|ρ=R2=0(R1<s<R2).Near an eigenvalueλi, the Green function G(ρ,s;λ) can be presented in the form (see, e.g., [22, 23])
(19)G(ρ,s;λ)=∑i=1∞vi(ρ)vi(s)λ-λi=vi(ρ)vi(s)λ-λi+G1(ρ,s;λ),
where the series converges absolutely and uniformly for all ρ,s∈[R1,R2] and for all λ from any compact set, which does not contain λi, i=1,2,…. The function G1(ρ,s;λ) is regular in the neighborhood of λi, while λn,vn(ρ) are the above-mentioned systems of eigenvalues and normalized in L2(R1,R2;ρ) eigenfunctions, see [24]. For big numbers n it is true that the eigenvalue λn=O(n2) (see [22]).Using the same technique as in [13], we obtain an integral representation for the solution u(s) of (11) for s∈[R1,R2]:
(20)u(s)=α∫R1R2G(ρ,s)ρu3(ρ)dρ+f(s),
where f(s)=R2u′(R2-0)G(R2,s)-R1u′(R1+0)G(R1,s).Assume thatf∈C[R1,R2] and λ≠λn. From the properties of the Green function, it follows that the kernel N(ρ,s;λ) is a continuous function in the square R1⩽ρ, s⩽R2.It follows from the results of the paper [13] that the nonlinear operator
(21)S(u)∶=∫R1R2N(ρ,s;λ)u3(ρ)dρ+f(s),
where N(ρ,s;λ)∶=αG(ρ,s;λ)ρ, is completely continuous on each bounded subset of C[R1,R2].Consider the cubic equation∥N∥r3+∥f∥=r, where
(22)∥N∥=max⁡s∈[R1,R2]∫R1R2|N(ρ,s)|dρ,∥f∥=max⁡s∈[R1,R2]|f(s)|.It can be shown (see [16]) that the following assertions hold, which will be needed below.Proposition 1.
Under the condition(23)0<∥f(λ)∥<2313∥N(λ)∥.
The equation ∥N∥r3+∥f∥=r has two nonnegative solutions
(24)r*=-213∥N∥cos⁡(13arccos(332∥f∥∥N∥)-2π3),r*=-213∥N∥cos⁡(13arccos(332∥f∥∥N∥)+2π3),
such that r*<r*.Proposition 2.
Under the condition∥f∥⩽(2/3)(1/3∥N∥), (20) has at least one solution such that ∥u∥⩽r*.Using Propositions1 and 2, it can be shown (see [16]) that the following theorems hold.Theorem 3.
Ifα⩽A2, where
(25)A=231∥f∥3∥N1∥,∥N1∥=max⁡∫R1R2|ρG(ρ,s)|dρ,
then (20) has a unique continuous solution: u∈C[R1,R2] such that ∥u∥⩽r*.Theorem 4.
Let the kernelN and the right-hand side f of integral equation (20) be continuous functions of the parameter  λ∈Λ0, N(ρ,s;λ)⊂C(Λ0×[R1,R2]×[R1,R2]), f(s,λ)⊂C(Λ0×[R1,R2]), on a certain interval Λ0 of the real axis. Let also
(26)0<∥f(λ)∥<2313∥N(λ)∥.
Then, the solutionu(ρ;λ) to (20), for λ∈Λ0, exists, is unique, and continuously depends on the parameter λ, u(ρ;λ)∈C(Λ0×[R1,R2]).We should note that the statement of the theorem holds for sufficiently smallα.Approximate solutionsun of integral equation (20), written in the form u=S(u), can be found using the iterative process un+1=S(un), n=0,1,2,…, with u0=0. The iterative procedure together with the convergence theorem is given in [13]. These iterated solutions can be applied to determine approximated eigenvalues of problem P.It should be noted that there is a simpler and more efficient approach to determine approximate eigenvalues and eigenfunctions which is based on the Cauchy problem method [20, 25] (shooting method) with initial data taken, for example, at the point ρ=R1. Of course, when performing the calculations, we need to specify the constant C1; in fact, this means specifying the eigenfunction at the point ρ=R1. As was noticed previously, the eigenvalues in the nonlinear problem depend on the initial conditions. The iterative algorithm described above yields an approximate solution but is too cumbersome for implementation. The main purpose of the theoretical apparatus used here is to rigorously prove the existence of eigenvalues. The numerical results presented below were obtained by the Cauchy problem method.


## 6. Dispersion Equation

Settings=R1+0 and s=R2-0 in (20), using transmission conditions (14), and taking into account solutions (12) and (13), we obtain
(27)C1I1(k1R1)=α∫R1R2G(ρ,R1)ρu3(ρ)dρ+C2R2k3K1′(k3R2)G(R2,R1)-C1R1k1I1′(k1R1)G(R1,R1),C2K1(k3R2)=α∫R1R2G(ρ,R2)ρu3(ρ)dρ+C2R2k3K1′(k3R2)G(R2,R2)-C1R1k1I1′(k1R1)G(R1,R2).Expressing the constantC2 from the first and second equations of the system and equating the resulting expressions, we obtain a dispersion equation of the form
(28)C1g(λ)=αF(λ),
where
(29)g(λ)=[I1(k1R1)+R1k1I1′(k1R1)G(R1,R1)]×[K1(k3R2)-R2k3K1′(k3R2)G(R2,R2)]+R1R2k1k3I1′(k1R1)K1′(k3R2)×G(R1,R2)G(R2,R1),F(λ)=[K1(k3R2)-R2k3K1′(k3R2)G(R2,R2)]×∫R1R2G(ρ,R1;λ)ρu3(ρ)dρ+R2k3K1′(k3R2)G(R2,R1)×∫R1R2G(ρ,R2;λ)ρu3(ρ)dρ.


## 7. Existence of Solutions of the Dispersion Equation

The zeros of the functionΦ(λ)≡C1g(λ)-αF(λ) are the values of λ for which problem P has a nontrivial solution; that is, if λ=λ~ is such that Φ(λ~)=0, then the eigenvalues of the problem are determined by the equation λ~=γ2.Let us analyze the existence of solutions of the dispersion equation for the linear problem. From (28), with α=0, we obtain g(λ)=0. This equation can be rewritten as
(30)-f1f2=f1G(R2,R2)+f2G(R1,R1)+G(R1,R1)G(R2,R2)-G(R1,R2)G(R2,R1),
where
(31)f1=I1(k1R1)R1k1I1′(k1R1)>0,f2=-K1(k3R2)R2k3K1′(k3R2)>0,
and we took into account that I1(k1R1)·K1(k3R2)≠0.Using (19) one finds
(32)G(R1,R1)=vi2(R1)λ-λi+G1(R1,R1),G(R2,R2)=vi2(R2)λ-λi+G1(R2,R2).
As series in (19) converges absolutely and uniformly it is possible to reorder terms of the series. So, after some algebra one can find that
(33)G(R1,R1)G(R2,R2)-G(R1,R2)G(R2,R1)=∑i=1∞‍∑j=1∞‍vi2(R1)vj2(R2)-vi(R1)vi(R2)vj(R1)vj(R2)(λ-λi)(λ-λi)=∑i=1∞‍∑j>i∞‍[vi(R1)vj(R2)-vi(R2)vj(R1)]2(λ-λi)(λ-λi)=∑i=1∞‍∑j>i∞‍[vi(R1)vj(R2)-vi(R2)vj(R1)]2λj-λimmmmmmmm×(-1λ-λi+1λ-λj).It is supposed that the eigenvaluesλi are ordered in the way that
(34)λi>λi+1,foranyi
and so, it is true that λj-λi<0 if j>i.Using derived formula (30) can be rewritten in the following form
(35)-f1f2-f1G1(R2,R2)-f2G1(R1,R1)=ai2f1+bi2f2λ-λi-∑i=1∞‍∑j>i∞‍cij2(-1λ-λi+1λ-λj),
where ai2=vi2(R2)≠0, bi2=vi2(R1)≠0, and
(36)cij2=[vi(R1)vj(R2)-vi(R2)vj(R1)]2(λi-λj)-1⩾0.
The left-hand side of (35) is bounded and continuous in the vicinity of the pole λi.Let us consider the behavior of the right-hand side of (35) in the vicinity of the same pole.Asf1(λi)+f2(λi)≠0 one can easily see that
(37)lim⁡λ→λi±0(-ai2f1+bi2f2λ-λi-∑i=1∞‍∑j>i∞‍cij2(-1λ-λi+1λ-λj))=±∞.It follows from this consideration that at least one root of the equationg(λ)=0 lies between adjacent poles (or, what the same is, adjacent eigenvalues of problem (16)) λi and λi+1.The fact that the termvi(Rp) in the expression G(Rp,Rq;λ), where p,q∈{1,2}, does not vanish can be easily proved using the theorem of the uniqueness of a solution to the Cauchy problem [13].Now, we can show that the solutions of the equation(38)Φ(λ)=0
exist.Indeed, let there be integersk⩾1 and l⩾0, such that
(39)max⁡(ε1,ε3)<λl<λl+1<⋯<λl+k-1<λl+k<ε~2,
where ε~2=min⁡ρ∈[R1,R2]ε2(ρ); λl+i, i=0,k¯ are eigenvalues of problem (16). Above we rearranged eigenvalues in the order of magnitude.Choose sufficiently small numbersδi>0 such that the Green function G(ρ,s;λ) exists and is continuous on the union
(40)Γ∶=⋃i=1k‍Γi
of the intervals
(41)Γi∶=[λi-1+δi-1,λi-δi],i=1,k¯.In addition, we assume thatδi are chosen such that the inequalities
(42)g(λi-1+δi-1)g(λi-δi)<0,i=i,k¯
are valid.In other words, we choseδi in such a way which move away from the poles of the Green function but preserve solutions of the linear problem (for α=0).Hence,F(λ) is bounded. Moreover, the product αF(λ) can be made sufficiently small by choosing a suitable α. Consider the dispersion equation Φ(λ)=0. The function g(λ) is continuous and changes its sign when λ varies from λi-1+δ~i-1 to λi-δ~i. Since F(λ) is bounded for λ varying from λi-1+δ~i-1 to λi-δ~i, by choosing a suitable α, we can always make the equation Φ(λ)=0 have at least k roots λ~i, such that λ~i∈(λi-1+δ~i-1,λi-δ~i), i=1,k¯, where δ~i are chosen on δi.The main result of the present work is the following theorem.Theorem 5.
Let the numbersε1, ε3, and ε~2=min⁡ρ∈[R1,R2]⁡ε2(ρ) satisfy the condition ε~2>max⁡(ε1,ε3)>0 and let there be integers k⩾1 and l⩾0 such that
(43)max⁡(ε1,ε3)<λl<λl+1<⋯<λl+k-1<λl+k<ε~2,
where λi are the eigenvalues of boundary value problem (16). Then, there is a number α0>0 such that, for any α⩽α0, problem P has at least k eigenvalues γi such that γi∈(λi-1+δi-1,λi-δi), i=1,k¯.Proof.
The Green function exists for allγ∈Γ. Clearly,
(44)A(γ)=23∥f(γ)∥3∥N1(γ)∥
is a continuous function of γ∈Γ. Let A1=min⁡γ∈Γ⁡A(γ) and let α<A12. According to Theorems 3, and 4, there is a unique solution u=u(γ) to equation (13) for any γ∈Γ. This solution is a continuous function, and ∥u∥⩽r*=r*(γ). Let r00=max⁡γ∈Γ⁡r*(γ). Estimating F(λ), we obtain |F(λ)|⩽Cr003, where C is a constant.
The functiong(γ) is continuous, and the equation g(γ)=0 has at least one root γ~i inside the interval Γi, λi-1+δi-1<γ~i<λi-δi. Let M1=min⁡l⩽i⩽l+k-1⁡|g(λi+δi)|, M2=min⁡l+1⩽i⩽l+k⁡|g(λi-δi)|. Then, M~=min⁡{M1,M2} is positive and independent of α.
Ifα⩽M~/Cr003, then
(45)(C1g(λi-1+δi-1)-αF(λi-1+δi-1))×(C1g(λi-δi)-αF(λi-δi))<0.
The functionC1g(λ)-αF(λ) is continuous; hence, the equation C1g(λ)-αF(λ)=0 has a root γ^i inside Γi; that is, λi+δi<γi<λi+1-δi+1. We can choose α0=min⁡{A12,M~/Cr003}.Theorem5 implies that, under the above-formulated conditions, there exist symmetric TE-polarized waves propagating without attenuation in nonhomogeneous cylindrical dielectric circular waveguides filled with a nonmagnetic isotropic Kerr nonlinear medium. This result generalizes the well-known assertion for both homogeneous dielectric waveguides (see [3]) and nonhomogeneous dielectric circular waveguides (see [4]) filled with a linear medium (with α=0).


## 8. Existence of Propagation Constants

Consider the Cauchy problem for(46)u′′=-ρ-1u′+ρ-2u-k2(ρ)u-αu3
with the initial conditions
(47)u(R1)∶=u(R1+0),u′(R1)∶=u′(R1+0),
where the values u(R1+0) and u′(R1+0) are determined from (12) and have the following forms:
(48)u(R1+0)=C1I1(k1R1),u′(R1+0)=C1k1(I0(k1R1)-I0(k1R1)k1R1),
where C1 is supposed to be known.It is clear from the foregoing that the calculation necessitates specified constantC1 or C2. The constant can be specified on any of the waveguide boundaries. Note that, in contrast to problems for linear media, the propagation constants (eigenvalues) in problems for nonlinear media depend on the constant at one of the boundaries.We use the classic results of the theory of ordinary differential equations on the existence and uniqueness of the solution to the Cauchy problem and on the continuous dependence of the solution to the Cauchy problem on the parameter [24].Letmax⁡(ε1,ε3)<γ*<γ*<∞ and let bγ<∞ be a certain constant. We determine the set
(49)Πγ:={(ρ,u;γ):|ρ-R1|⩽ρ~,|u-u(R1)|⩽bγ,γ∈[γ*,γ*]},
and number Mγ is such that Mγ⩾max⁡Πγ⁡|P|, where P is the vector of right-hand side of (46) if it is rewritten in normal form.The following propositions are valid (see [26] for the proof).Proposition 6.
Solutionu(ρ;γ) to Cauchy problem (46)-(47) exists for all ρ∈[R1,R2], where R2⩽min⁡(ρ~,bγ/Mγ), is unique, continuously depends on γ for all γ∈[γ*,γ*], and is continuously differentiable with respect to ρ.Remark 7.
Since we consider the caseR1>0, it is evident that Proposition 6 has a nonlocal character; that is, we can always choose R2 and γ*,γ* such that solution u(ρ;γ) exists and is continuous for ρ∈[R1,R2] and for all γ∈[γ*,γ*].Let, forγ=γ~, the equality
(50)u(R2-0;γ~)=u(R2+0;γ~)
be fulfilled. If the equality
(51)u′(R2-0;γ~)=u′(R2+0;γ~)
is also fulfilled, then γ~ is an eigenvalue of problem P.Inherently, formulae (50) and (51) are transmission conditions (14) at the point ρ=R2.In order to determineγ~, suppose that (50) is fulfilled and consider the function
(52)F(R2;γ)∶=u′(R2-0;γ)-u′(R2+0;γ).From transmission conditions (14) on the boundary ρ=R2 and solution (13), we obtain
(53)u(R2-0)=C2K1(k3R2),u′(R2-0)=-C2k3(K0(k3R2)+K1(k3R2)k3R2).Using this system, we can rewriteF(R2;γ) in the following form:
(54)F(R2;γ)=u′(R2-0;γ)+k3u(R2-0;γ)K1(k3R2)(K0(k3R2)+K1(k3R2)k3R2).It follows from this formula that the valueF(R2;γ) is expressed in terms of values of the solution to the Cauchy problem only.One can show [26, 27] that the following statement is valid.Proposition 8.
Assume that the conditions of Proposition6 are fulfilled and the segment [γ_,γ¯]⊂[γ*,γ*] is such that F(R2;γ_)·F(R2;γ¯)<0. Then, there is, at least, one propagation constant (one eigenvalue) γ~∈(γ_,γ¯) of problem P.Remark 9.
The conditionF(R2;γ_)·F(R2;γ¯)<0 is sufficient for the existence of propagation constant γ~∈(γ_,γ¯) of problem P.It follows from general results of the theory of ordinary differential equations that calculated eigenvalues of the problem are stable with respect to initial values.


## 9. Determination of Approximate Propagation Constants

The method under consideration makes it possible to plot the dependence of (normalized) propagation constantγ on the (normalized) radius ΔR=R2-R1 of the waveguide. Curves γ=γ(ω) (or γ=γ(f)), where ω=2πf is the circular frequency, are dispersion curves in such problems. If curve γ(ω) depends on the incident field amplitude (a circumstance we deal with in the considered problem), these curves are referred to as energy dispersion curves [6]. Since normalized variables are employed, the plot of the dependence γ=γ(ΔR) is called a dispersion curve (or an energy dispersion curve).Consider Cauchy problem (46)-(47). Let 0<ΔR*<ΔR*<∞, max⁡(ε1,ε3)<γ*<γ*<∞ be certain numbers. Assume that,
(55)ΔR∈[ΔR*,ΔR*],γ∈[γ*,γ*].We split the segments[ΔR*,ΔR*] and [γ*,γ*] into n and m sections, respectively. Then, we have {ΔRi,γi}, i=0,n¯ and j=0,m¯. Here, ΔR0=ΔR*, ΔRn=ΔR*, γ0=γ*, and γm=γ*. Then, for each pair of indices (i,j), we have pair of initial values (uij(R1),uij′(R1)), where uij(R1)∶=u(R1+0;γj) and uij(R1)∶=u′(R1+0;γj). In spite of the fact that the initial values are independent of ΔRi, the double indexing is suitable for our analysis.Now, we can formulate the Cauchy problem for (46) with the initial condition uij(R1),uij′(R1). The solution of this problem yields the values
(56)u(ΔRi;γj)∶=uij(R2-0)u′(ΔRi;γj)∶=uij′(R2-0).Suppose thatuij(R2-0)=uij(R2+0) and construct the function
(57)F(ΔRi;γj)∶=uij′(R2-0)-uij′(R2+0).Assume that, for givenΔRi, there exist γj and γj+1 such that
(58)F(ΔRi;γj)·F(ΔRi;γj+1)<0.Hence, there existsγ~j∈(γj,γj+1) that is an eigenvalue of the considered problem of wave propagation, and this eigenvalue is associated with the layer’s thickness ΔRi. Value γ~j can be found within an arbitrary accuracy by means of, for example, the dichotomy method.On the basis of the dichotomy method, we develop a technique for the determination of an approximate value of the propagation constant.Let us specify errorη>0 for the determination of the value of propagation constant γ~. Let the segment (γ_1,γ¯1) be such that
(59)F(ΔRi,γ_1)·F(ΔRi,γ¯1)<0.The sought eigenvalue isγ~∈(γ_1,γ¯1) and the approximate eigenvalue is γ~1∈(γ_1,γ¯1).Let us determine the center of the segmentγ1=0.5(γ_1+γ¯1) and calculate F(ΔRi;γ1). We check the following conditions.(1)
If|F(ΔRi;γ1)|<η, then γ1 is the sought approximate eigenvalue.(2)
IfF(ΔRi;γ_1)·F(ΔRi;γ1)<0, then γ~∈(γ_1,γ1). In this case, we set γ_2∶=γ_1 and  γ¯2∶=γ1; hence, γ~2∈(γ_2,γ¯2).(3)
IfF(ΔRi;γ1)·F(ΔRi;γ¯1)<0, then γ~∈(γ1,γ¯1). In this case, we set γ_2∶=γ1 and  γ¯2∶=γ¯1, hence; γ~2∈(γ_2,γ¯2).We repeat dichotomic divisionp times to obtain the approximate value γ~p∈(γ_p,γ¯p). Evidently, |γ¯p-γp|=2-p|γ¯1-γ1|. Let us choose number p such that 2-p|γ¯1-γ1|<η. Then, for example, the center of the segment (γ_p,γ¯p), that is, γ~p=0.5(γ¯p+γ_p), can be considered as approximate value γ~p of γ~ where γ~.One can show [26] that the following statement is valid.Proposition 10.
Assume that the conditions of Proposition8 are fulfilled and that, by means of the dichotomy method, sequence {γ~p} of approximate values of propagation constant γ~ is obtained; then, lim⁡p→∞⁡γ~p=γ~, where γ~ is an exact eigenvalue.


## 10. Numerical Results

Numerical results are found with the help of the method suggested in the previous section. Numerical experiments were carried out for two types of nonlinearities: the Kerr law and nonlinearity with saturation. For the inhomogeneityɛ(ρ) of the waveguide, the following functions, that specify the permittivities ɛ=ɛ(ρ)+αu2 and ɛ=ɛ(ρ)+αu2/(1+βu2) in the layer R1<ρ<R2, are used.(1)
Forɛ=ɛ(ρ)+αu2,
(a)
ɛ(ρ)=ε2+1/ρ;(b)
ɛ(ρ)=ε2 (homogeneous case);(c)
ɛ(ρ)=ε2+ρ.(2)
Forɛ=ɛ(ρ)+αu2/(1+βu2),
(a)
ɛ(ρ)=ε2+1/ρ;(b)
ɛ(ρ)=ε2 (homogeneous case);(c)
ɛ(ρ)=ε2+ρ,whereε2 is a positive real constant.In figures below eigenvalues or dependenceγ(ΔR) (see Figures 2–4 and 7–9) and eigenfunctions u(ρ) (see Figures 5, 6, 10, and 11) are shown.Figure 2
Case 1(a). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.524 (red dot γ2≈6.370) and nonlinear eigenvalues γ≈2.573 and γ≈2.962 (green and black dots, resp., γ2≈6.620 and γ2≈8.770).Figure 3
Case 1(b). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.461 (red dot γ2≈6.056) and nonlinear eigenvalues γ≈2.491 and γ≈2.999 (green and black dots, resp., γ2≈6.200 and γ2≈8.990).Figure 4
Case 1(c). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.619 (red dot γ2≈6.859) and nonlinear eigenvalues γ≈2.712 and γ≈2.939 (green and black dots, resp., γ2≈7.355 and γ2≈8.638).Figure 5
Eigenfunctions for Case 1(a). The inner radiusR1=2 and the outer radius R2=3. Red line (linear problem) γ≈2.524, green line (nonlinear problem) γ≈2.573, and black line (nonlinear problem) γ≈2.962.Figure 6
Eigenfunctions for Case 1(c). The inner radiusR1=2 and the outer radius R2=2.75. Red line (linear problem) γ≈2.619, green line (nonlinear problem) γ≈2.712, and black line (nonlinear problem) γ≈2.939.Figure 7
Case 2(a). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.524 (red dot γ2≈6.370) and nonlinear eigenvalues γ≈2.571, γ≈3.142, and γ≈3.806 (green, blue, and black dots, resp., γ2≈6.610, γ2≈9.872, and γ2≈14.486).Figure 8
Case 2(b). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.461 (red dot γ2≈6.056) and nonlinear eigenvalues γ≈2.490, γ≈3.219, and γ≈3.736 (green, blue, and black dots, resp., γ2≈6.200, γ2≈10.362, and γ2≈13.958).Figure 9
Case 2(c). The inner radiusR1=2. Marked eigenvalues: linear eigenvalue γ≈2.619 (red dot γ2≈6.859) and nonlinear eigenvalues γ≈2.701, γ≈3.131, and γ≈3.779 (green, blue, and black dots, resp., γ2≈7.295, γ2≈9.803, and γ2≈14.281).Figure 10
Eigenfunctions for Case 2(a). The inner radiusR1=2 and the outer radius R2=3. Red line (linear problem) γ≈2.524, green line (nonlinear problem) γ≈2.571, blue line (nonlinear problem) γ≈3.142, and black line (nonlinear problem) γ≈3.806.Figure 11
Eigenfunctions for Case 2(c). The inner radiusR1=2 and the outer radius R2=2.75. Red line (linear problem) γ≈2.619, green line (nonlinear problem) γ≈2.701, blue line (nonlinear problem) γ≈3.131, and black line (nonlinear problem) γ≈3.779.In fact, Cases 1(b) and 2(b) are homogeneous; however, we consider them as inhomogeneous with constant inhomogeneity (it prevents us from mentioning each time the difference).The following values of parameters are used for calculations given in Figures2–11: ε1=4, ε2=9, ε3=1, C1=1, and R1=2.For the Kerr nonlinearity (Figures2–6), α=0.01.For the nonlinearity with saturation (Figures7–11), α=0.01, and β=0.001.Dispersion curves for Cases 1(a), 1(b), and 1(c) are shown in Figures2–4, respectively; dispersion curves for Cases 2(a), 2(b), and 2(c) are shown in Figures 7–9, respectively.In Figures2–4 and Figures 7–9, the vertical axis corresponds to γ2>4 and the horizontal axis corresponds to the thickness 0<ΔR=R2-R1<6 of the outer cylindrical shell (see Figure 1).In Figures5, 6, 10, and 11, the vertical axis corresponds to the value of the plotted function u(ρ) and the horizontal axis corresponds to the value ρ∈[0,5].Red lines in Figures2–4 and 7–9 correspond to the linear inhomogeneous case and blue curves correspond to the nonlinear inhomogeneous case (in the same figure, the inhomogeneity is the same for both types of curves; red and blue curves differ from each other only by the nonlinear term).Vertical dashed line in Figures2, 3, 7, and 8 corresponds to R2-R1=1; in Figures 4 and 9  vertical dashed line corresponds to R2-R1=0.75. Eigenvalues of the problem are points of intersections of this dashed line with the dispersion curves; first three points in Figures 2–4 are marked (the smallest value, red dot, corresponds to the linear case and the others two (green and black dots) correspond to the Kerr case); first four points in Figures 7–9 are marked (the smallest value, red dot, corresponds to the linear case and the other three (green, blue, and black dots) correspond to the case with saturated nonlinearity).In Figures5 and 6, eigenfunctions for Cases 1(a) and 1(c), respectively, are plotted. The color of a curve corresponds to the color of marked eigenvalue in Figures 2 and 4. The larger an eigenvalue is, the higher maximum corresponding eigenfunction has.In Figures10 and 11, eigenfunctions for the Cases 2(a) and 2(c), respectively, are plotted. The color of a curve corresponds to the color of marked eigenvalue in Figures 7 and 9. The larger an eigenvalue is, the higher maximum corresponding eigenfunction has.If we compare red (linear case) and blue (nonlinear case) dispersion curves in Figures2–4, we can notice that in the linear cases there is only finite number of the eigenvalues (this fact is well known) and for the Kerr case is we can suppose infinite number of the eigenvalues (as it is for the case of plane waveguide with, the Kerr nonlinearity; see, for example, [16, 28]). In each nonlinear case, there are eigenvalues that are close to corresponding linear case (these eigenvalues can be determined with the help of perturbation theory; existence of these very eigenvalues is proved analytically in this paper). At the same time, it is easy to see, from Figures 2–4 and 7–9, that in each nonlinear case there are eigenvalues that cannot be determined from perturbation theory (we call them “purely nonlinear” eigenvalues).Calculation time usually is about 2–5 hours with a modern laptop (it depends on the grid density).


## 11. Conclusion and Discussion

As it is obvious, the analytical method allows proving existence of eigenvalues in the nonlinear problem which are close to eigenvalues in the corresponding linear problem. However, as it is seen from the plots above, in the nonlinear problem, there are new eigenvalues that cannot be determined as a perturbation of eigenvalues in the linear problem. These new eigenvalues correspond to the new (“purely” nonlinear) propagation regime. At the same time, proposed numerical method allows finding all eigenvalues. Of course this method is efficient in the case of discrete eigenvalues.In Figures2–4, groups of three eigenvalues are marked. The red dot (smallest value) is an eigenvalue of the linear problem (with α=0). The green dot is an eigenvalue of the nonlinear problem. When one passes to the limit α→0, then the green dot tends to the red dot. The black dot does not tend to an eigenvalue of the linear problem, when α→0; this eigenvalue is “purely” nonlinear and corresponds to a new type of nonlinear waves.The same conclusion can be made for the eigenvalues shown in Figures7–9. In these figures, groups of four eigenvalues are marked. Green, blue, and black dots correspond to the nonlinear problem. The green dot is a perturbation of the solution of the linear problem (red dot) and tends to the red dot, when α→0. Blue and black dots are “purely” nonlinear eigenvalues; they also correspond to a new type of nonlinear waves.Whether these mathematically predicted “purely” nonlinear waves really exist, it is a hypothesis that can be proved or disproved in the physical experiment.The numerical method considered in the paper has the following advantages.(+) The method is easy to be implemented, because the Cauchy problem for a system of ordinary differential equations can be solved by means of the standard tools of any mathematical software package.(+) The method works substantially (“substantially” here means that the implementation of the numerical method based on auxiliary Cauchy problem is simpler and this method calculates faster approximately by a factor of ten.) faster than the numerical method based on the use of integral dispersion equations.(+) The method allows determining of approximate eigenvalues within an arbitrary given accuracy.And the following is its disadvantage.(−) Whenγ~ is one of the eigenvalues of the problem, then the total derivative of function F(ΔRi;γ) will respect to γ must not vanish at γ=γ~.This method is also applied for plane-multilayered structures; see the paper “electromagnetic wave propagation in nonlinear-layered waveguide structures: computational approach to determine propagation constants” of Valovik in [29].



---

*Source: 231498-2014-01-22.xml*