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

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


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][2][3][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][10][11][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].

Statement of the Problem
Let us consider a three-dimensional space R 3 with cylindrical coordinate system .The space is filled with an isotropic medium of dielectric permittivity  3  0 = const, where  0 is 2 Advances in Numerical Analysis 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  := {(, , ) : 0 ⩽  <  1 , 0 ⩽  < 2} ∪ {(, , ) :  1 ⩽  ⩽  2 , 0 ⩽  < 2} (1) and is placed in R 3 with its generating line parallel to the axis .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 radii  1 and  2 ; that is, the waveguide is two-layered (see Figure 1):  1 is the radius of the inner cylinder and  2 −  1 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].
The geometry of the problem is shown in Figure 1.The waveguide is unlimitedly continued in  direction.
Let E and H be complex amplitudes of an electromagnetic field.The complex amplitudes E, H must satisfy Maxwell's equations [7] rot have continuous tangential field components on the media interface  =  1 ,  =  2 and obey the radiation condition at infinity; that is, the electromagnetic field decays exponentially as  → ∞ in the region  >  2 ;  is the circular frequency.
The permittivity in the entire space has the form  = ε 0 , where and e  is the orthonormal vector in the  direction; (⋅, ⋅) is the Euclidean scalar product;

Transmission Conditions and Transmission Problem
Transmission conditions for the functions  and   result from the continuity conditions for the tangential field components and have the form where is the jump of the limit values of the function at a point .
Let us formulate the transmission eigenvalue problem (problem ) 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  1 ̸ = 0 (or  2 ̸ = 0), there is a nonzero function (; ) that is defined by formulas ( 12) and (13) for  <  1 and  >  2 , respectively, and solves equation (11) for  1 <  <  2 ; moreover, the function (; ) thus defined satisfies conditions (8) and transmission conditions (14).
The quantities solving problem  are called eigenvalues, and the corresponding functions (; ) 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.
Let  :=  2 .Consider the boundary value problem Let   , V  () 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  ̸ =   , the boundary value problem has only the trivial solution.This means that, for  ̸ =   , there uniquely exists the Green function (, ; ) of the boundary value problem Near an eigenvalue   , the Green function (, ; ) can be presented in the form (see, e.g., [22,23]) where the series converges absolutely and uniformly for all ,  ∈ [ 1 ,  2 ] and for all  from any compact set, which does not contain   ,  = 1, 2, . ... The function  1 (, ; ) is regular in the neighborhood of   , while   , V  () are the above-mentioned systems of eigenvalues and normalized in  2 ( 1 ,  2 ; ) eigenfunctions, see [24].For big numbers  it is true that the eigenvalue   = ( 2 ) (see [22]).
Using the same technique as in [13], we obtain an integral representation for the solution () of (11) for  ∈ [ 1 ,  2 ]: where Assume that  ∈ [ 1 ,  2 ] and  ̸ =   .From the properties of the Green function, it follows that the kernel (, ; ) is a continuous function in the square  1 ⩽ ,  ⩽  2 .
It follows from the results of the paper [13] that the nonlinear operator where (, ; ) := (, ; ), is completely continuous on each bounded subset of [ 1 ,  2 ].
Consider the cubic equation ‖N‖ 3 + ‖‖ = , where It can be shown (see [16]) that the following assertions hold, which will be needed below.
Using Propositions 1 and 2, it can be shown (see [16]) that the following theorems hold.
then (20) has a unique continuous solution: Theorem 4. Let the kernel  and the right-hand side f of integral equation (20) be continuous functions of the parameter ), on a certain interval Λ 0 of the real axis.Let also Then, the solution (; ) to (20), for  ∈ Λ 0 , exists, is unique, and continuously depends on the parameter , We should note that the statement of the theorem holds for sufficiently small .
Approximate solutions   of integral equation ( 20), written in the form  = (), can be found using the iterative process  +1 = (  ),  = 0, 1, 2, . .., with  0 = 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 .
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  =  1 .Of course, when performing the calculations, we need to specify the constant  1 ; in fact, this means specifying the eigenfunction at the point  =  1 .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.

Dispersion Equation
Setting  =  1 + 0 and  =  2 − 0 in (20), using transmission conditions (14), and taking into account solutions ( 12) and ( 13), we obtain Expressing the constant  2 from the first and second equations of the system and equating the resulting expressions, we obtain a dispersion equation of the form where

Existence of Solutions of the Dispersion Equation
The zeros of the function Φ() ≡  1 () − () are the values of  for which problem  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 () = 0.This equation can be rewritten as where and we took into account that  1 ( 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 It is supposed that the eigenvalues   are ordered in the way that and so, it is true that Using derived formula (30) can be rewritten in the following form where The left-hand side of ( 35) is bounded and continuous in the vicinity of the pole   .Let us consider the behavior of the right-hand side of (35) in the vicinity of the same pole.
As  1 (  ) +  2 (  ) ̸ = 0 one can easily see that lim It follows from this consideration that at least one root of the equation () = 0 lies between adjacent poles (or, what the same is, adjacent eigenvalues of problem ( 16))   and  +1 .
The fact that the term V  (  ) in the expression (  ,   ; ), where ,  ∈ {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 exist.Indeed, let there be integers  ⩾ 1 and  ⩾ 0, such that where ε2 = min ∈[ 1 , 2 ]  2 ();  + ,  = 0,  are eigenvalues of problem (16).Above we rearranged eigenvalues in the order of magnitude.
Choose sufficiently small numbers   > 0 such that the Green function (, ; ) exists and is continuous on the union of the intervals In addition, we assume that   are chosen such that the inequalities are valid.
In other words, we chose   in such a way which move away from the poles of the Green function but preserve solutions of the linear problem (for  = 0).
The main result of the present work is the following theorem.
Theorem 5 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).

Existence of Propagation Constants
Consider the Cauchy problem for with the initial conditions where the values ( 1 +0) and   ( 1 +0) are determined from (12) and have the following forms: where  1 is supposed to be known.It is clear from the foregoing that the calculation necessitates specified constant  1 or  2 .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].
Let √max( 1 ,  3 ) <  * <  * < ∞ and let   < ∞ be a certain constant.We determine the set and number   is such that   ⩾ max Π  ||, where  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).Let, for  = γ, the equality is also fulfilled, then γ is an eigenvalue of problem .
One can show [26,27] that the following statement is valid.
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.

Determination of Approximate Propagation Constants
The method under consideration makes it possible to plot the dependence of (normalized) propagation constant  on the (normalized) radius Δ =  2 −  1 of the waveguide.Curves  = () (or  = ()), where  = 2 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  = (Δ) is called a dispersion curve (or an energy dispersion curve).Consider Cauchy problem (46)-(47 We split the segments [Δ * , Δ * ] and [ * ,  * ] into  and  sections, respectively.Then, we have {Δ  ,   },  = 0,  and  = 0, .Here, Δ 0 = Δ * , Δ  = Δ * ,  0 =  * , and   =  * .Then, for each pair of indices (, ), we have pair of initial values (  ( 1 ),    ( 1 )), where   ( 1 ) := ( 1 + 0;   ) and   ( 1 ) :=   ( 1 + 0;   ).In spite of the fact that the initial values are independent of Δ  , the double indexing is suitable for our analysis.Now, we can formulate the Cauchy problem for (46) with the initial condition   ( 1 ),    ( 1 ).The solution of this problem yields the values Suppose that   ( 2 − 0) =   ( 2 + 0) and construct the function Assume that, for given Δ  , there exist   and  +1 such that Hence, there exists γ ∈ (  ,  +1 ) that is an eigenvalue of the considered problem of wave propagation, and this eigenvalue is associated with the layer's thickness Δ  .Value γ 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 The sought eigenvalue is γ ∈ ( ,  1 ) and the approximate eigenvalue is γ1 ∈ ( 1 ,  1 ).Let us determine the center of the segment  1 = 0.5( 1 +  1 ) and calculate (Δ  ;  1 ).We check the following conditions.
One can show [26] that the following statement is valid.
Proposition 10.Assume that the conditions of Proposition 8 are fulfilled and that, by means of the dichotomy method, sequence {γ  } of approximate values of propagation constant γ is obtained; then, lim  → ∞ γ = γ, where γ is an exact eigenvalue.

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  = () +  2 and  = () +  2 /(1 +  2 ) in the layer  1 <  <  2 , are used.
( 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).
In Figures 5,6,10,and 11, the vertical axis corresponds to the value of the plotted function () and the horizontal axis corresponds to the value  ∈ [0, 5].
Red lines in Figures 2-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).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 Figures 5 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 Figures 10 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 Figures 2-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).

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 Figures 2-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 Figures 7-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 (Δ  ; ) 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].

Figure 1 :
Figure 1: Geometry of the problem.