# Nonlinear Effects of Electromagnetic TM Wave Propagation in Anisotropic Layer with Kerr Nonlinearity

Authors: Yu G. Smirnov; D. V. Valovik
Journal: Advances in Mathematical Physics
 (2012)

---

## Abstract

The problem of electromagnetic TM wave propagation through a layer with Kerr nonlinearity is considered. The layer is located between two half-spaces with constant permittivities. This electromagnetic problem is reduced to the nonlinear boundary eigenvalue problem for ordinary differential equations. It is necessary to find eigenvalues of the problem (propagation constants of an electromagnetic wave). The dispersion equation (DE) for the eigenvalues is derived. The DE is applied to nonlinear metamaterial as well. Comparison with a linear case is also made. In the nonlinear problem there are new eigenvalues and new eigenwaves. Numerical results are presented.


---

## Body

## 1. Introduction

Problems of electromagnetic wave propagation in nonlinear waveguide structures are intensively investigated during several decades. First known studies about nonlinear optics' problems are given in the monographs [1, 2]. Propagation of electromagnetic wave in a layer and a circle cylindrical waveguide are among such problems. Phenomena of electromagnetic wave propagation in nonlinear media have original importance and also find a lot of applications, for example, in plasma physics, microelectronics, optics, and laser technology. There are a lot of different nonlinear phenomena in media when an electromagnetic wave propagates, such as self-focusing, defocusing, and self-channeling [1–5].Investigation of nonlinear phenomena leads us to solve nonlinear differential equations. In some cases it is necessary to solve nonlinear boundary eigenvalue problems (NBEPs), which rarely can be solved analytically. One of the important nonlinear phenomenan is the case when the permittivity of the sample depends on electric field intensity. And one of the simplest nonlinearities is a Kerr nonlinearity [4, 6, 7]. When we speak about NBEPs we mean that differential equations and boundary conditions nonlinearly depend on the spectral parameter and also the differential equations nonlinearly depend on the sought-for functions. These facts do not allow to apply well-known methods of spectral problems' investigation.Here we consider electromagnetic TM wave propagation in a layer with Kerr nonlinearity. Perhaps, the papers [6, 7] were the first studies where some problems of electromagnetic wave propagation are considered in a strong electromagnetic statement. Propagation of polarized electromagnetic waves in a layer and in a circle cylindrical waveguide with Kerr nonlinearity is considered in this paper. When one says that the permittivity ɛ is described by Kerr law this means that (for an isotropic material) ɛ=ɛconst+α|E|2, ɛconst is the constant part of the permittivity ɛ; α is the nonlinearity coefficient; |E|2=Ex2+Ey2+Ez2, where E=(Ex,Ey,Ez) is an electric field. Below we consider an anisotropic case (some results were presented in [8]). The first approximation for eigenvalues of the problem is presented in [9].Problems of electromagnetic wave propagation in a linear layer (with constant permittivity) and in a linear circle cylindrical waveguide were deeply studied many years ago, see, for example, [10]. Such problems are formulated as boundary eigenvalue problems for ordinary differential equations. Indeed, the main interest in this problem is the value of the spectral parameter (eigenvalues) which corresponds to the propagating wave. If an eigenvalue is known it is easy to solve differential equations numerically. Otherwise numerical methods cannot be successfully applied. However, in nonlinear cases it is often paid more attention to solve the differential equations (see, e.g., [11–13]). Though the first problem is to find eigenvalues therefore to find Des, from the mathematical standpoint the DE is an equation with respect to the spectral parameter. Analysis of this equation allows us to make conclusions about problem's solvability, eigenvalues' localization, and so forth. In most cases the equations of the problem cannot be integrated in an explicit form. Of course, if one has the explicit solutions of the differential equations it is easy to derive the DE. Therefore, when the equations cannot be integrated things do not go to a DE. However, in some cases the DE can be found in an explicit form and it is not necessary to have explicit solutions of differential equations.Let us discuss in detail the case of Kerr nonlinearity. The work [4] contains a wide range of details of third-order nonlinear electromagnetic TE and TM guided waves. Problems of surface wave propagation along the interface between two semi-infinite linear or/and nonlinear media were studied completely (see the results in [4]). At the same time we should notice that problems of wave propagation in a nonlinear layer that is located between two semi-infinite linear or/and nonlinear media are much more difficult than (and cannot be reduced to) the problems where surface waves are considered only at the interface between two semi-infinite linear or/and nonlinear media. Propagation of TE waves was more studied. The work [14] is devoted to the problem of electromagnetic wave propagation in a nonlinear dielectric layer with absorption and the case of Kerr nonlinearity is considered separately. One of the most interesting works about propagation of TE waves in a layered structure with Kerr nonlinearity is the paper [15]. Also the reader can see the work [16], where a layer with Kerr nonlinearity without absorption is considered.The case of TM wave propagation in a nonlinear medium is more complicated. This is due to the fact that two components of the electric field make the analysis much harder [17].In the work [18] a linear dielectric layer is considered. The layer is located between two half-spaces. The half-spaces are filled by nonlinear medium with Kerr nonlinearity. This problem for TE waves is solved analytically [19, 20]. For the TM case in [18] obtained DE is an algebraic equation. It should be noticed that in [18] authors simplify the problem. Earlier in [21] the DE is obtained with other simplifying assumption (authors take into account only one component Ex of the electric field). Later in [22] it is proved that the dominating nonlinear contribution in the permittivity is proportional to the transversal component Ez. In the works [11] propagation of TM waves in a nonlinear half-space with Kerr nonlinearity is considered. Formal solutions of differential equations in quadratures are obtained. In the paper [11] DEs are presented for isotropic and anisotropic media in a half-space with nonlinear permittivity. The DEs are rational functions with respect to the value of field's components at the interface. Authors found the first integral of the system of differential equations (so called a conservation law). This is also very interesting work to study, another way to simplify the problem pointed out in [23].In the case of TE wave you can see the papers [24–26]. Propagation of TM wave in terms of the magnetic component is studied in [12, 13]. The paper [21] is devoted to the question (from physical standpoint) why it is possible to take into account only one component of the electric field in the expression for permittivity in the case of TM waves in a nonlinear layer. The results are compared with the case of TE waves.The most important results about TM wave propagation in a layer with Kerr nonlinearity (system of differential equations, first integral) and a circle cylindrical waveguide (system of differential equations) were obtained in [6, 7]. In some papers (e.g., [12]) polarized wave propagation in a layer with arbitrary nonlinearity is considered. However, DEs were not obtained and no results about solvability of the boundary eigenvalue problem were obtained as well. The problem of TM wave propagation in a layer with Kerr nonlinearity is solved at first for a thin layer and then for a layer of arbitrary thickness [27–29]. Theorems of existence and localization of eigenvalues are proved in [30, 31]. Some numerical results are shown in [8, 9].In this paper the DE is an equation with additional conditions. Only for linear media (when permittivity is a constant) in a layer or in a circle cylindrical waveguide the DEs are sufficiently simple (but even for these cases the DEs are transcendental equations). For a nonlinear layer the DE is quite complicated nonlinear integral equation, where the integrand is defined by implicit algebraic function. It should be stressed that in spite of the fact that the DE is complicated it can be rather easily solved numerically.This DE allows to study both nonlinear materials and nonlinear metamaterials. It should be noticed that in this paper materials with nonlinear permittivity and constant positive permeability are studied. But it is not difficult to take into account the sign of the permeability.Problems of propagation of TE and TM waves in a nonlinear circle cylindrical waveguide are also close to the problem considered here. These problems are more complicated in comparison with corresponding problems in nonlinear layers. And even in the case of Kerr nonlinearity the results are not so complete as in the case in layers [30, 32, 33].


## 2. Statement of the Problem

Let us consider electromagnetic wave propagation through a homogeneous anisotropic nonmagnetic dielectric layer. The layer is located between two half-spaces:x<0 and x>h in Cartesian coordinate system Oxyz. The half-spaces are filled with isotropic nonmagnetic media without any sources and characterized by permittivities ɛ1≥ɛ0 and ɛ3≥ɛ0, respectively, where ɛ0 is the permittivity of free space. Assume that everywhere μ=μ0 is the permeability of free space (see Figure 1).Figure 1
The geometry of the problem.It should be noticed that conditionsɛ1≥ɛ0, ɛ3≥ɛ0 are not necessary. They are not used for derivation of DEs, but they are useful for DEs' solvability analysis.The electromagnetic field depends on time harmonically [6]:
(2.1)E~(x,y,z,t)=E+(x,y,z)cos⁡ωt+E-(x,y,z)sin⁡ωt,H~(x,y,z,t)=H+(x,y,z)cos⁡ωt+H-(x,y,z)sin⁡ωt,
where ω is the circular frequency; E~, E+, E-, H~, H+, H- are real functions. Everywhere below the time multipliers are omitted.Form complex amplitudes of the electromagnetic field(2.2)E=E++iE-,H=H++iH-,
where E=(Ex,Ey,Ez)T, H=(Hx,Hy,Hz)T, and (·)T denotes the operation of transposition, and each component of the fields is a function of three spatial variables.Electromagnetic field(E,H) satisfies the Maxwell equations
(2.3)rot⁡H=-iωɛE,rot⁡E=iωμH,
the continuity condition for the tangential field components on the media interfaces x=0, x=h, and the radiation condition at infinity: the electromagnetic field exponentially decays as |x|→∞ in the domains x<0 and x>h.The permittivity inside the layer is described by the diagonal tensor(2.4)ɛ^=(ɛxx000ɛyy000ɛzz),
where
(2.5)ɛxx=ɛ2x+b|Ex|2+a|Ez|2,ɛzz=ɛ2z+a|Ex|2+b|Ez|2
and a,b,ɛ2>max(ɛ1,ɛ3) are positive constants (below the solutions are sought under more general conditions). It does not matter what a form ɛyy has. Since ɛyy is not contained in the equations below for the TM case, it should be noticed that ɛ^ describes tensor Kerr nonlinearity. When a=b we obtain scalar Kerr nonlinearity. Moreover, chosen nonlinearity satisfies the condition ∂ɛxx/∂Ez2=∂ɛzz/∂Ex2. This equation is satisfied by almost every known nonlinear Kerr mechanism, such as electronic distortion, molecular orientation, electrostriction, and Kerr nonlinearities described within the uniaxial approximation mentioned in the paper [11]. The case when ɛ2x=ɛ2z is studied in [8]. Pay heed to the fact that the problem considered here is not studied in [31].The solutions to the Maxwell equations are sought in the entire space.


## 3.   TM Waves

Let us consider TM waves(3.1)E=(Ex,0,Ez)T,H=(0,Hy,0)T,
and Ex, Ez, Hy are functions of three spatial variables. It is easy to show that the components of the fields do not depend on y. Waves propagating along medium interface z depend on z harmonically. This means that the fields components have the form
(3.2)Ex=Ex(x)eiγz,Ez=Ez(x)eiγz,Hy=Hy(x)eiγz,
where γ is the spectral parameter of the problem.So we obtain from system (2.3) [6]
(3.3)iγEx(x)-Ez′(x)=iωμHy(x),Hy′(x)=-iωɛzzEz(x),iγHy(x)=iωɛxxEx(x),
where (·)′≡d/dx.The following equation can be easily derived from the previous system:(3.4)Hy(x)=1iωμ(iγEx(x)-Ez′(x)).Differentiating (3.4) and using the second and the third equations of system (3.3) we obtain
(3.5)γ(iEx(x))′-Ez′′(x)=ω2μɛzzEz(x),γ2(iEx(x))-γEz′(x)=ω2μɛxx(iEx(x)).Let us denotek02:=ω2μ0ɛ0 and perform the normalization according to the formulas x~=k0x, d/dx=k0(d/dx~), γ~=γ/k0, ɛ~j=ɛj/ɛ0 (j=1,2,3), a~=a/ɛ0, b~=b/ɛ0. Denoting by Z(x~):=Ez, X(x~):=iEx and omitting the tilde symbol, from system (3.5) we obtain(3.6)-Z′′+γX′=ɛzzZ,-Z′+γX=γ-1ɛxxX.It is necessary to find eigenvaluesγ of the problem that correspond to surface waves propagating along boundaries of the layer 0<x<h. We seek the real values of the spectral parameter γ such that real solutions X(x) and Z(x) to system (3.6) exist. Indeed, in this case |E|2 does not depend on z. Since E=(Ex(x)eiγz,0,Ez(x)eiγz)=eiγz(Ex,0,Ez), therefore, |E|=|eiγz|·|Ex|2+|Ez|2. It is known that |eiγz|=1 as Imγ=0. Let γ=γ′+iγ′′. Then, we obtain |eiγz|=|eiγ′z|·|e-γ′′z|=e-γ′′z. If γ′′≠0, then e-γ′′z is a function on z. In this case the components Ex and Ez depend on z, but it contradicts to the choice of Ex(x) and Ez(x). So we can consider only real values of γ.We consider that(3.7)ɛ={ɛ1,x<0,ɛ^,0<x<h,ɛ3,x>h.Also we assume thatmax(ɛ1,ɛ3)<γ2<min(ɛ2x,ɛ2z). This two-sided inequality naturally appears for the problem in a layer with a constant permittivity tensor.FunctionsX, Z are supposed to be sufficiently smooth due to physical nature of the problem
(3.8)X(x)∈C(-∞,0]∩C[0,h]∩C[h,+∞)∩C1(-∞,0]∩C1[0,h]∩C1[h,+∞),Z(x)∈C(-∞,+∞)∩C1(-∞,0]∩C1[0,h]∩C1[h,+∞)∩C2(-∞,0)∩C2(0,h)∩C2(h,+∞).It is clear that system (3.6) is an autonomous one. System (3.6) can be rewritten in a normal form. This system in the normal form can be considered as a dynamical system with analytical with respect to X and Z right-hand sides. Of course, in the domain where these right-hand sides are analytical with respect to X and Z, it is well known (see, e.g., [34]) that the solutions X and Z of such a system are analytical functions with respect to the independent variable as well. This is an important fact for DEs' derivation.We consider thatγ2>max(ɛ1,ɛ3).


## 4. Differential Equations of the Problem

In the domainx<0 we have ɛ=ɛ1. From system (3.6) we obtain X′=γZ, Z′=γ-1(γ2-ɛ1)X. In accordance with the radiation condition we obtain
(4.1)X(x)=Aexγ2-ɛ1,Z(x)=-Aγ-1γ2-ɛ1exγ2-ɛ1.We assume thatγ2-ɛ1>0; otherwise it will be impossible to satisfy the radiation condition.In the domainx>h we have ɛ=ɛ3. From system (3.6) we obtain X′=γZ, Z′=γ-1(γ2-ɛ3)X. In accordance with the radiation condition we obtain(4.2)X(x)=Be-(x-h)γ2-ɛ3,Z(x)=-Bγ-1γ2-ɛ3e-(x-h)γ2-ɛ3.Here for the same reason as above we consider thatγ2-ɛ3>0.ConstantsA and B in (4.1) and (4.2) are defined by transmission conditions and initial conditions.Inside the layer0<x<h system (3.6) takes the form(4.3)-d2Zdx2+γdXdx=[ɛ2z+aX2+bZ2]Z,-dZdx+γX=γ-1[ɛ2x+bX2+aZ2]X.Differentiating the second equation and substituting its right-hand side instead of left-hand side into the first equation we can rewritten system (4.3) in the following form:(4.4)dXdx=2aγɛ2x-γ2+bX2+aZ2ɛ2x+3bX2+aZ2X2Z+γɛ2z+aX2+bZ2ɛ2x+3bX2+aZ2Z,dZdx=-γ-1[ɛ2x-γ2+bX2+aZ2]X.Now system (4.4) is written in a normal form. If the right-hand sides are analytic functions with respect to X and Z, then the solutions are analytic functions with respect to its independent variable.Dividing the first equation in system (4.4) to the second one we obtain the ordinary differential equation
(4.5)-(ɛ2x+3bX2+aZ2)dXdZ=2aXZ+γ2ɛ2z+aX2+bZ2ɛ2x-γ2+bX2+aZ2ZX.Equation (4.5) can be transformed into a total differential equation.Its solution (first integral of system (4.4)) can be easily found and be written in the following form:
(4.6)C=b2X6+2abX4Z2+a2X2Z4+12(4ɛ2x-3γ2)bX4+(2ɛ2x-γ2)aX2Z2+12γ2bZ4+γ2(ɛ2x-γ2)X2+(ɛ2x-γ2)2X2+γ2ɛ2zZ2,
where C is a constant of integration.


## 5.   Transmission Conditions and the Transmission Problem

Tangential components of an electromagnetic field are known to be continuous at media interfaces. In this case the tangential components areHy and Ez. Hence, we obtain
(5.1)Hy(h+0)=Hy(h-0),Hy(0-0)=Hy(0+0),Ez(h+0)=Ez(h-0),Ez(0-0)=Ez(0+0).From the continuity conditions for the tangential components of the fieldsE and H and using (3.4) we obtain
(5.2)γX(h)-Z′(h)=Hy(h),γX(0)-Z′(0)=Hy(0),Z(h)=Ez(h+0)=Ez(h),Z(0)=Ez(0-0)=Ez(0),
where Hy(h):=i(μ/ɛ0)Hy(h+0), Hy(0):=i(μ/ɛ0)Hy(0-0).The constantEz(h):=Ez(h+0) is supposed to be known (initial condition). Let us denote X0:=X(0), Xh:=X(h), Z0:=Z(0), and Zh:=Z(h). So we obtain that A=(γ/γ2-ɛ1)Z0, B=(γ/γ2-ɛ3)Zh.Then from conditions (5.2) we obtain
(5.3)Hy(h)=-Zhɛ3γ2-ɛ3,Hy(0)=Z0ɛ1γ2-ɛ1.
In accordance with (3.6), (3.7) inside the layer
(5.4)-Z′(x)+γX(x)=γ-1(ɛ2x+bX2(x)+aZ2(x))X(x).
Then for x=h, using (5.2), we obtain from (5.4)
(5.5)γ-1[ɛ2x+bXh2+aZh2]Xh=Hy(h).
From (5.5) we obtain the equation with respect to Xh:
(5.6)Xh3+b-1(ɛ2x+aZh2)Xh-b-1γHy(h)=0.Under taken assumptions (in regard toɛ2 and a) the value a-1(ɛ2+aZh2)>0. Hence, this equation has at least one real root, which is considered (the root can be find explicitly by using Cardanus-Ferrari formula [35].Using first integral (4.6) at x=h, we find the value ChX:=C|x=h from the equation
(5.7)ChX=b2Xh6+2abXh4Zh2+a2Xh2Zh42-1(4ɛ2x-3γ2)bXh4+(2ɛ2x-γ2)aXh2Zh2+2-1γ2bZh4+γ2(ɛ2x-γ2)Xh2+(ɛ2x-γ2)2Xh2+γ2ɛ2zZh2.In order to find the valuesX0 and Z0 it is necessary to solve the following system (this system is obtained using formula (5.4) at x=0 and the first integral at the same point):
(5.8)γɛ1Z0=γ2-ɛ1[ɛ2x+bX02+aZ02]X0,ChX=b2X06+2abX04Z02+a2X02Z04+2-1(4ɛ2x-3γ2)bX04+(2ɛ2x-γ2)aX02Z02+2-1γ2bZ04+γ2(ɛ2x-γ2)X02+(ɛ2x-γ2)2X02+γ2ɛ2zZ02.It is easy to see from the second equation of this system that the valuesX0 and Z0 can have arbitrary signs. At the same time from the first equation of this system we can see that X0 and Z0 must be positive or negative simultaneously.Normal components of electromagnetic field are known to be discontinues at media interfaces. And it is the discontinuity of the first kind. In this case the normal component isEx. It is also known that the value ɛEx is continuous at media interfaces. From the above and from the continuity of the tangential component Ez it follows that the transmission conditions for the functions ɛX and Z are (5.9)[ɛX]x=0=0,[ɛX]x=h=0,[Z]x=0=0,[Z]x=h=0,
where [f]x=x0=limx→x0-0f(x)-limx→x0+0f(x) denotes a jump of the function f at the interface.We also suppose that functionsX(x) and Z(x) satisfy the condition
(5.10)X(x)=O(1|x|),Z(x)=O(1|x|)as|x|⟶∞.


## 6. Dispersion Equation

Introduce the new variables(6.1)τ(x)=ɛ2x+bX2(x)+aZ2(x)γ2,η(x)=γX(x)Z(x)τ(x).Using new variables rewrite system (4.4),
(6.2)dτdx=2τη(γ2τ-ɛ2x)bη2+aγ2τ2×(bη2+aγ2τ2)(bɛ2z-aγ2τ(τ-1))+(aη2+bγ2τ2)(γ2τ-ɛ2x)γ2τ(bη2+aγ2τ2)+2bη2(γ2τ-ɛ2x),dηdx=τ-1τη2+ɛ2z+(γ2τ-ɛ2x)aη2+bγ2τ2bη2+aγ2τ2,
and (4.6) (6.3)γ2τ-ɛ2xbη2+aγ2τ2[η2((γ2τ-ɛ2x)2+ɛ2x(ɛ2x-γ2))+γ4ɛ2zτ2]+(γ2τ-ɛ2x)22(bη2+aγ2τ2)2×[(4ɛ2x-3γ2)bη4+2(2ɛ2x-γ2)aγ2τ2η2+γ6bτ4]=C,
where constant C is equal to the constant C in (4.6).In order to obtain the DE for the propagation constants it is necessary to find the valuesη(0), η(h).It is clear thatη(0)=γ(X(0)/Z(0))τ(0), η(h)=γ(X(h)/Z(h))τ(h). Taking into account that γ2X(x)τ(x)=ɛX(x) and using formulas (5.2), (5.3), it is easy to obtain that
(6.4)η(0)=ɛ1γ2-ɛ1>0,η(h)=-ɛ3γ2-ɛ3<0.It is easy to see that the right-hand side of the second equation of system (VI) is strictly positive. This means that the functionη(x) monotonically increases on interval (0,h). Taking into account (6.4) we obtain that the function η(x) cannot be differentiable on the entire interval (0,h). This means that the function η(x) has a break point. Let x*∈(0,h) be the break point. From (6.3) it is obvious that x* is such that τ*=τ(x*) is a root of the equation Chτ+3(τ*)2-2(τ*)3-2τ0(2-τ*)τ*=0. In addition η(x*-0)→+∞ and η(x*+0)→-∞.It is natural to suppose that the functionη(x) on interval (0,h) has several break points x0,x1,…,xN. The properties of function η(x) imply(6.5)η(xi-0)=+∞,η(xi+0)=-∞,wherei=0,N¯.Let(6.6)1w:=τ-1τη2+ɛ2z+(γ2τ-ɛ2x)aη2+bγ2τ2bη2+aγ2τ2,
where w=w(η); τ=τ(η) is expressed from (5.4).Taking into account our hypothesis we will seek the solutions on each interval[0,x0),(x0,x1),…,(xN,h]:(6.7)-∫η(x)η(x0)‍wdη=x+c0,∫η(xi)η(x)‍wdη=x+ci,∫η(xN)η(x)‍wdη=x+cN,
where 0≤x≤x0, xi≤x≤xi+1, and xN≤x≤h, respectively, and i=0,N-1¯.Substitutingx=0, x=xi+1, and x=xN into equations in (6.7) (into the first, the second, and the third, resp.,) and taking into account (6.5), we find constants c1,c2,…,cN+1:
(6.8)c0=-∫η(0)+∞‍wdη,ci+1=∫-∞+∞‍wdη-xi+1,cN+1=∫-∞η(h)‍wdη-h,
where i=0,N-1¯.Using (6.8) we can rewrite (6.7) in the following form:
(6.9)∫η(x)η(x0)‍wdη=-x+∫η(0)+∞‍wdη,∫η(xi)η(x)‍wdη=x+∫-∞+∞‍wdη-xi+1,∫η(xN)η(x)‍wdη=x+∫-∞η(h)‍wdη-h,
where 0≤x≤x0, xi≤x≤xi+1, xN≤x≤h, respectively, and i=0,N-1¯.Introduce the notationT:=∫-∞+∞‍wdη. It follows from formula (6.9) that 0<xi+1-xi=T<h, where i=0,N-1¯. This implies the convergence of the improper integral (it will be proved in other way below). Now consider x in (6.9) such that all the integrals on the left side vanish (i.e., x=x0, x=xi, and x=xN), and sum all equations in (6.9). We obtain
(6.10)0=-x0+∫η(0)+∞‍wdη+x0+T-x1+⋯+xN-1+T-xN+xN+∫-∞η(h)‍wdη-h.Finally we obtain(6.11)-∫η(h)η(0)‍wdη+(N+1)T=h,
where η(0), η(h) are defined by formulas (6.4).Expression (6.11) is the DE, which holds for any finite h. Let γ be a solution of DE (6.11) and an eigenvalue of the problem. Then, there are eigenfunctions X and Z, which correspond to the eigenvalue γ. The eigenfunction Z has N+1 zeros on the interval (0,h).Notice that improper integrals in DE (6.11) converge. Indeed, function τ=τ(η) is bounded as η→∞ since τ=γ-2(ɛ2x+bX2+aZ2), and X, Z are bounded.Then(6.12)|w|≤1αη2+β,
where α>0, β>0 are constants. It is obvious that improper integral ∫-∞+∞‍(dη/(αη2+β)) converges. Convergence of the improper integrals in (6.11) in inner points results from the requirement that the right-hand side of the second equation of system (VI) is positive.The first equation of system (VI) jointly with the first integral can be integrated in hyperelliptic functions. The solution is expressed in implicit form by means of hyperelliptic integrals. This is the simple example of Abelian integrals. The inversion of these integrals is hyperelliptic functions and they are solutions of system (VI). Hyperelliptic functions are Abelian functions, which are meromorphic and periodic functions. Since functionη is expressed algebraically through τ, therefore, η is a meromorphic periodic function. This means that the break point x* is a pole of function η.


## 7. Generalized Dispersion Equation

Here we derive the generalized DE, which holds for any real valuesɛ2. In addition the sign of the right-hand side of the second equation in system (VI) and condition max(ɛ1,ɛ3)<γ2<ɛ2 are not taken into account. These conditions appear in the case of a linear layer and are used for derivation of DE (6.11). Though on the nonlinear case it is not necessary to limit the value γ2 from the right side, at the same time it is clear that γ is limited from the left side, since this limit appears from the solutions in the half-spaces.Now we assume thatγ satisfies the following two-sided inequality: (7.1)max⁡(ɛ1,ɛ3)<γ2<+∞.Using first integral (6.3) it is possible to integrate formally any of the equations of system (VI). As earlier we integrate the second equation, we cannot obtain the solution on the entire interval (0,h), since function η(x) can have break points, which belong to (0,h). It is known that function η(x) has break points only of the second kind (η is an analytical function).Assume that functionη(x) on interval (0,h) has N+1 break points x0,x1,…,xN.It should be noticed that(7.2)η(xi-0)=±∞,η(xi+0)=±∞,
where i=0,N¯, and signs ± are independent and unknown.Taking into account the previous, solutions are sought on each interval[0,x0),(x0,x1),…,(xN,h]:
(7.3)-∫η(x)η(x0-0)‍wdη=x+c0,∫η(xi+0)η(x)‍wdη=x+ci+1,∫η(xN+0)η(x)‍wdη=x+cN+1,
where 0≤x≤x0, xi≤x≤xi+1, and xN≤x≤h, respectively, and i=0,N-1¯.From (7.3), substituting x=0, x=xi+1, and x=xN into the first, the second, and the third equations in (7.3), respectively, we find required constants c1,c2,…,cN+1:
(7.4)c0=-∫η(0)η(x0-0)‍wdη,ci+1=∫η(xi+0)η(xi+1-0)‍wdη-xi+1,cN+1=∫η(xN+0)η(h)‍wdη-h,
where i=0,N-1¯.Using (7.4), (7.3) take the form
(7.5)∫η(x)η(x0-0)‍wdη=-x+∫η(0)η(x0-0)‍wdη,∫η(xi+0)η(x)‍wdη=x+∫η(xi+0)η(xi+1-0)‍wdη-xi+1,∫η(xN+0)η(x)‍wdη=x+∫η(xN+0)η(h)‍wdη-h,
where 0≤x≤x0, xi≤x≤xi+1, and xN≤x≤h, respectively, and i=0,N-1¯.From formulas (7.5) we obtain that(7.6)xi+1-xi=∫η(xi+0)η(xi+1-0)‍wdη,i=0,N-1¯.Expressions0<xi+1-xi<h<∞ imply that under the assumption about the break point existence the integral on the right side converges and ∫η(xi+0)η(xi+1-0)‍wdη>0. In the same way, from the first and the last equations of (7.5) we obtain that x0=∫η(0)η(x0-0)‍wdη and 0<x0<h then (7.7)0<∫η(0)η(x0-0)‍wdη<h<∞,
and h-xN=∫η(xN+0)η(h)‍wdη and 0<h-xN<h then
(7.8)0<∫η(0)η(x0-0)‍wdη<h<∞.These considerations yield that the functionw(η) has no nonintegrable singularities for η∈(-∞,∞). And also this proves that the assumption about finite number break points is true.Now, settingx=x0, x=xi, and x=xN into the first, the second, and the third equations in (7.5), respectively, we have that all the integrals on the left sides vanish. We add all the equations in (7.5) to obtain
(7.9)0=-x0+∫η(0)η(x0-0)‍wdη+x0+∫η(x0+0)η(x1-0)‍wdη-x1+⋯+xN-1+∫η(xN-1+0)η(xN-0)‍wdη-xN+xN+∫η(xN+0)η(h)‍wdη-h.From (7.9) we obtain(7.10)∫η(0)η(x0-0)‍wdη+∫η(xN+0)η(h)‍wdη+∑i=0N-1‍∫η(xi+0)η(xi+1-0)‍wdη=h.It follows from formulas (7.6) that
(7.11)η(xi+0)=±∞,η(xi-0)=∓∞,wherei=0,N¯,
and it is necessary to choose the infinities of different signs.Thus we obtain that(7.12)∫η(x0+0)η(x1-0)‍wdη=⋯=∫η(xN-1+0)η(xN-0)‍wdη=:T′.
Hence x1-x0=⋯=xN-xN-1.Now we can rewrite (7.10) in the following form: (7.13)∫η(0)η(x0-0)‍wdη+∫η(xN+0)η(h)‍fdη+NT′=h.LetT≡∫-∞+∞‍wdη; then we finally obtain
(7.14)-∫η(h)η(0)‍wdη±(N+1)T=h,
where η(0), η(h) are defined by formulas (6.4).Expression (7.14) is the DE, which holds for any finite h. Let γ be a solution of DE (7.14) and an eigenvalue of the problem. Then, there are eigenfunctions X and Z, which correspond to the eigenvalue γ. The eigenfunction Z has N+1 zeros on the interval (0,h). It should be noticed that for every number N+1 it is necessary to solve two DEs: for N+1 and for -(N+1).Note 1.
If there is a certain valueγ*2, such that some of the integrals in DEs (6.11) or (7.14) diverge at certain inner points this simply means that the value γ*2 is not a solution of chosen DE and the value γ*2 is not an eigenvalue of the problem.Note 2.
It is necessary to emphasize that this boundary eigenvalue (transmission) problem essentially depends on the initial conditionZh. The transmission problem for a linear layer does not depend on the initial condition. If the nonlinearity function is a specific one, then in some cases it will be possible to normalize the Maxwell equations in such a way that the transmission problem does not depend on initial condition Zh explicitly (it is possible e.g., for Kerr nonlinearity in a layer and in a circle cylindrical waveguide). Once more we stress the fact that the opportunity of such normalization is an exceptional case. What is more, in spite of the fact that in certain cases this normalization is possible it does not mean that the normalized transmission problem is independent of the initial condition. In this case one of the problem's parameter depends on the initial condition.


## 8. Passage to the Limit in the Generalized Dispersion Equation

In this section we assume thatɛ2x=ɛ2z=ɛ2 and b=a. Now consider the passage to the limit as a→0. The value a=0 corresponds to the case of a linear medium in the layer. Two cases are possible:(a)ɛ2>0,(b)ɛ2<0 (metamaterial case).Let us examine case (a). The DE for a linear case is well known [10] and has the form
(8.1)tg(hɛ2-γ2)=ɛ2ɛ2-γ2(ɛ1γ2-ɛ3+ɛ3γ2-ɛ1)ɛ1ɛ3(ɛ2-γ2)-ɛ22γ2-ɛ3γ2-ɛ1.Let(8.2)f=τγ2τ2+η2(τ-1),f1=ɛ2ɛ2-γ21ɛ22/(ɛ2-γ2)+η2.Using passage to the limit asa→0 we obtain the function f1 from the function f. We seek bounded solutions X(x) and Z(x). This implies that the denominator of the function f1 cannot vanish. What is more, the function f as a→0 tends to the function f1 uniformly on x∈[0,h]. It is possible to pass to the limit under integral sign as a→0 in (7.14) using results of classical analysis
(8.3)h=-ɛ2ɛ2-γ2∫η(h)η(0)‍1ɛ22/(ɛ2-γ2)+η2dη+ɛ2ɛ2-γ2(N+1)∫-∞+∞‍1ɛ22/(ɛ2-γ2)+η2dη,
where η(0), η(h) are defined by formulas (6.4).The integrals in (8.3) are calculated analytically. Calculating these integrals we obtain
(8.4)hɛ2-γ2=arctgɛ2ɛ2-γ2(ɛ1γ2-ɛ3+ɛ3γ2-ɛ1)ɛ1ɛ3(ɛ2-γ2)-ɛ22γ2-ɛ3γ2-ɛ1+(N+1)π.Expression (8.4) can be easily transformed into expression (8.1).Let us examine (b) case. We haveɛ2<0 (metamaterial) and the DE for the linear case has the form [31]
(8.5)e2hγ2-ɛ2=ɛ1γ2-ɛ2-ɛ2γ2-ɛ1ɛ1γ2-ɛ2+ɛ2γ2-ɛ1ɛ3γ2-ɛ2-ɛ2γ2-ɛ3ɛ3γ2-ɛ2+ɛ2γ2-ɛ3,
where γ2-ɛ1>0, γ2-ɛ2>0, and γ2-ɛ3>0.In the same way as above, passing to the limit in the functionf as a→0 we obtain f2=|ɛ2|/(γ2-ɛ2)(1/η2-ɛ22/(γ2-ɛ2)). Passing to the limit in equation (7.14) as a→0 and integrating the function f2, after simple calculations, we obtain formula(8.5)The results here show that it is possible to pass to the limit asa→0. DE in (7.14) for the nonlinear case turns into (8.1) or (8.5) for the linear case as a→0.


## 9. Numerical Results

The way of solution to the DE in (7.14) is the following: we choose the segment on γ then cut this segment into p pieces with nods γi,i=0,1,…,p. Then for each γi we can calculate all necessary values in (7.14). The integral in (7.14) can be calculated using any method of numerical integration. In order to calculate the value of the integrand at the point γi it is necessary to use first integral (6.3). For each value γi we calculate value hi, so we obtain the grid {γi,hi},i=0,1,…,p. Choosing reasonably dense grid on γ we can plot dependence γ on h, as it is done below.Dispersion curves (DC) calculated from (7.14), (8.1), and (8.5) are shown in Figures 2 and 4. Eigenmodes for eigenvalues indicated in Figures 3 and 5 are shown in Figures 3 and 5.Figure 2
Plot ofγ(h). The first few dispersion curves are shown. Solid curves for the nonlinear case (solutions of (7.14)); dashed curves for the linear case (solutions of (8.1)). The following parameters are used for both cases: ɛ1=1.44, ɛ2=9, ɛ3=1, and for the nonlinear case a=0.1, and Ez(h)=1. Dashed lines are described by formulas: h*=3.206 (thickness of the layer), γ=1.2 (lower bound for γ), and γ=3 (upper bound for γ in the case of linear medium in the layer).Eigenfunctions (fields) for the nonlinear problem are shown. Solid curves forX; dashed curves for Z. The same parameters as in Figure 2 are used. For (a), γ=2.994; for (b), γ=3.892: for (c), γ=8.657, and h=3.206 is used for all three cases. The eigenvalues are marked in Figure 2.
(a)
(b)
(c)Figure 4
Plot ofγ(h). The first few dispersion curves are shown. Solid curves for the nonlinear case (solutions of (7.14)); dashed curve for the linear case (solutions of (8.5)). The following parameters are used for both cases: ɛ1=1, ɛ2=-1.5, ɛ3=1, and for the nonlinear case a=5.2 and Ez(h)=1. Dashed lines are described by formulas: h*=1.71 (thickness of the layer), γ=1 (lower bound for γ) and γ=3 (upper bound for γ in the case of linear medium in the layer).Eigenfunctions (fields) for the nonlinear problem are shown. Solid curves forX; dashed curves for Z. The same parameters as in Figure 4 are used. For (a), γ=2.620, and a=0 (see (8.5)); for (b), γ=1.565; for (c), γ=3.481, and h=1.71 is used for all three cases. The eigenvalues are marked in Figure 4.
(a)
(b)
(c)As it is known and it is shown in Figure2, the line γ=3 is an asymptote for DCs in the linear case. It should be noticed that in the linear case there are no DCs in the region γ2≥ɛ2. It can be proved that function h≡h(γ) defined from equation (7.14) is continuous at the neighborhood γ2=ɛ2 when a≠0 (see Figure 2). This is the important distinction between linear and nonlinear cases.Further, it can be proved that functionh≡h(γ) defined from equation (7.14) when a≠0 has the following property: (9.1)lim⁡γ2→+∞h(γ)=0.In Figure2 for h=3.206 in the case of a linear layer there are 3 eigenvalues (black dots where the line h=3.206 intersects with DCs). These eigenvalues correspond to 3 eigenmodes. In the case of a nonlinear layer in Figure 2 are shown 5 eigenvalues (uncolored dots). These eigenvalues correspond to 5 eigenmodes. Taking into account the last paragraph's statement it is clear that in this case there is infinite number of eigenvalues.It is easy to see in Figure6 that the less nonlinearity coefficient a is the more stretched DCs in the nonlinear case. The maximum points of the curves h(γ) (in Figure 6 they are marked by asterisks) move to the right. The parts of the DCs that locate below the maximum points asymptotically tend to the DCs for the linear case as a→0.Figure 6
Plot ofγ(h) for the different values of a: 1 – a=1; 2 – a=0.1; 3 – a=0.01; 4 – a=0.001; 5 – a=0.0001; 6 – a=0 (linear case). The following parameters are used for both cases: ɛ1=4, ɛ2=9, ɛ3=1, and for the nonlinear case Ez(h)=1. Dashed lines are described by formulas: γ=2 (lower bound for γ), γ=3 (upper bound for γ in the case of linear medium in the layer). Curves (solid) 1–5 are solutions of (7.14), and curve 6 (dashed) is  solutions of (8.1).It should be noticed that in the case of Kerr nonlinearity in a layer and TE waves there are strong constraints on the valuea depending on the value ɛ2 (for details see [16]). It is natural to suppose that there are constraints of the kind in the case under consideration.As far as we know experiments to observe the new nonlinear eigenmodes were not carried out. So the question if the modes corresponding to the new eigenvalues exist (in an experiment) stays open! We should like to emphasize that it is interesting to observe purely nonlinear modes that do not arise in a linear limiting case. If to see Figure2 then points (b) and (c) correspond to these purely nonlinear modes.



---

*Source: 609765-2012-07-01.xml*