Guided TE-waves in a slab structure with lossless cubic nonlinear dielectric and magnetic material: parameter dependence and power flow with focus on metamaterials

The parameter dependence and power flow of guided TE-waves in a lossless cubic nonlinear, dielectric, magnetic planar three-layer structure is studied as follows. Using a travelling wave ansatz with stationary amplitude, Maxwell’s equations are transformed to a system of ordinary nonlinear differential equations.The solutions of the system are presented compactly (in terms of hyperbolic and elliptic functions).The nonnegative and bounded (“physical”) solutions are determined by using a phase diagram condition (PDC) that is applied to express the continuity (transmission) conditions at the interfaces leading to the dispersion relation (DR).Based on the PDC, the parameter dependence and stability of the solutions to the DR and corresponding power flow are studied numerically for permittivities and permeabilities that may be appropriate to describe metamaterial.


Introduction
The theory of electromagnetic wave propagation in nonlinear dielectric slab structure usually considers only nonmagnetic dielectrics [1][2][3]. Obviously, to be applied to metamaterials, the theory must be extended by including magnetic material. Assuming a linear (with respect to the magnetic field H) permeability this extension is straightforward (see Sects. 2 and 3), and most of relevant articles devoted to guided waves in metamaterial slab structures [4][5][6][7][8][9][10][11][12] are considering a nonlinear permittivity and linear permeability.
Apart from different particular assumptions for the material parameters in cladding and substrate (see Fig.  1), three of above articles consider a metamaterial film [4,6,7], thus making a certain contact with the present article in this respect.
In [4], Darmanyan et al. are investigating waves in a slab wave guide with negative index nonlinear film surrounded by "traditional" linear semi-infinite media. They are using particular solutions (Eqs. (1) in [4]) to obtain special dispersion relations and special expressions for the power flow. The conditions the problems' parameters must satisfy to use equations (2a)-(2d) are not presented. Thus, the scope of applications seems rather restricted.
Based on earlier publications, Boardman and Egan [6] are studying waves in a lossless, cubic nonlinear, double-negative film bounded by standard, lossless, linear, non-dispersive, non-magnetic dielectric cladding and substrate. Depending on different problems' parameters different Jacobi functions are used as solutions of the Helmholtz equation (see Eq. (1) in [6]), leading to numerical evaluation of the dispersion relations, field profiles, and power flow.
Compared with the present paper, there are two major differences. The first is the use of a phase diagram [13] to exclude nonphysical (nonreal and unbounded) solutions of the Helmholtz equation from the beginning. The second is to use the Weierstrass' elliptic functions instead of Jacobi elliptic functions as solutions of the Helmholtz equation. In principle, this difference is marginal. However, since we are interested in the parameter dependence of the solutions, variable problems' parameters (and hence variable propagation constants as solutions of the dispersion relation) imply variable modulus, and hence different Jacobi functions, and hence different representations of the dispersion relation. Thus, numerical evaluation becomes involved without a compact representation of the field intensities and of the dispersion relation.
Manenkov [7] is considering complex nonlinear permittivity and permeability µ, depending on the transverse coordinate in the film. For substrate and cladding complex linear and µ are assumed. Maxwell's equations are solved numerically. Since we use a different model (see Eqs. (1)-(3) below) it seems appropriate to do without a comparison here.
We note that none of references [4][5][6][7][8][9][10][11][12] deals with parameter dependence of solutions to the problem (stated in the next section). But parameter dependence is crucial for practical applications. If, e.g., the thickness of the waveguide is not suitably related to the material parameters, the solution of Maxwell's equations is singular in the central layer (see Fig. 6). Similarly, the same holds for the incident intensity J 0 and the material parameters.
The method for studying parameter dependence proposed below is based on the analytical properties of the solutions J(x, γ) (see Eq. (4)) combined with the dependence of the derivatives dJ dx 2 on J, represented as a phase diagram (see Fig. 2). It is essential for the straightforward (though involved) evaluation that function J f (x, γ; p j ), expressed in terms of Weierstrass function ℘(x; g 2 , g 3 ) (see [13]), analytically depends on the parameters p j . Hence, evaluation of the dispersion relation (14) is simple (using MATEMATICA with built-in function ℘). By selecting the nonnegative and bounded intensities J f by considering the corresponding phase diagrams (see Fig. 2), the dependence of solution γ of (14) on the parameters p j can be studied easily (see Fig. 3).
The paper is organized as follows. Section 2 describes the problem. Section 3 is devoted to the solutions of the basic nonlinear differential equations, to the phase diagram condition (PDC), to the dispersion relation together with its solvability, and, finally, to the parameter dependence of the propagation constant γ and of the total power flow P tot . By selecting particular parameters, the results of Section 3 are elucidated in Section 4. The paper concludes with comments in Section 5.

Statement of the problem
The structure we consider is shown in Figure 1. The propagation down the guide is in the z-direction. The guide is lying in the (z, y)-plane so that the problem is independent on y. We assume negligible small loss, isotropic, homogeneous, and magnetic material in the layers with local permittivities depending on the transverse electric field E y (x, z, t) as with ν , a ν , ν = s, f, c, real and constant, and permeabilities µ according to with real constants µ ν and µ 0 as the free space permeability. We note here, that the model described by equations (1) and (2) can be used for traditional material with , µ > 0. For metamaterial ( , µ < 0) it is a rather wide spread approach despite the problems connected with homogenization of the periodic (meta) structures as mentioned in the Introduction. Returning to this problem below, we disregard it for the present, and start with a tentative (TE)-solution (withẼ y , γ real) that transforms (see Eqs. (3) and (4) in [13]) Maxwell's equations subject to equations (1) and (2) to the system of ordinary nonlinear differential equations See equation (4) above.
The problem is to find physical (nonnegative and bounded) solutions J(x) to equation (4) that satisfy the continuity (transmission) conditions at the interface x = 0 and x = h together with the radiation conditions at infinity In particular, the problem is to find physical solutions J ν (x, γ) for metamaterial layers so that dispersion relations and power flow can be expressed in terms of J ν (x, γ) and (J ν (x, γ)), suitable to study the parameter dependence (and thus stability) of the solutions J ν (x, γ).

Solution
Apart from J(x, γ; p j ) = constant(p j denotes the fixed parameters of the problem), the solutions of equation (4) are [13]: See equations (6) and (7) above.
where J 0 , J (h) denote the intensities at x = 0 and x = h, respectively. The invariants g 2 , g 3 of Weierstrass' elliptic function ℘(x; g 2 , g 3 ) are given by with C f to be determined by the transmission conditions at x = 0 or at x = h (see below). Integration constants C s , C c are zero according to conditions (5). The prime in (7) denotes differentiation w.r.t. x for ℘(x; g 2 , g 3 ) and differentiation w.r.t. J for R f (J).
Continuity of E y and E y µ at the surfaces between the layers implies continuity of J and J µ (transmission conditions). Using equations (6)-(9), evaluation of J 0 = J s± (0, γ; p j ) = J f± (0, γ; p j ) and Rs(J0) Similarly, equations can be evaluated to give a restriction of the (unknown) intensity J (h) in equation (10), Due to the continuity of J at x = h, the continuity condition for ∂xJ µ at x = h can be written as Equation (14) represents the dispersion relation (DR) for a waveguide with permittivities and permeabilities given by (1) and (2), respectively. Though equation (13) and the dispersion relation are coupled by requirement it is useful to consider them separately. Since we are interested in physical solutions J, at least one root of equation (13) must be positive. If the discriminant D of equation (13) is nonnegative, and then (first case) only one positive root J (h) exists or (second case) two positive roots exist, respectively. Ambiguity of the latter case can be removed by the requirement J (h) (γ, p j ) = J f± (h, γ; p j ), since physical J f± (h, γ; p j ), existence provided, are defined uniquely according to equation (8). In this context we note that, in general, and if R f (J 0 ) > 0, but unbounded for vanishing denominator in equations (7)- (9). It can be seen from equations (7) or (8), J f± (x p , γ; p j ) = J 0 at the poles x p of ℘ and ℘ . The constraint R f (J 0 ) > 0, that is necessary for real J f± (x, γ; p j ) and ∂ x J f± (x, γ; p j ), can be satisfied only for particular propagation constants γ and particular parameters { ν , µ ν , a ν , J 0 }. Due to the transmission conditions, The qualitative behaviour of solutions J ν of equation (4) ((J ) 2 = R ν (J)) conveniently can be determined by considering phase diagrams {R ν (J), J} [14] as depicted in Figure 2. Since R ν (J) ≥ 0 must hold, with J varying monotonically until J = 0, it is clear that the zeros of R ν (J) are important. Some thought shows that ten and only ten phase diagrams must be taken into account (phase diagrams (a)-(j) for J f , phase diagram (k) for J s,c ) to characterize physical (real, nonnegative, bounded) solutions J ν . In order to find solutions γ = γ (h; p j ) of the dispersion relation (14), first, J f± must be physical (for x ∈ (0, h)) with R c (J ± (h, γ; p j )) > 0, and second, the dispersion relation must be solvable. These requirements conveniently can be studied by considering the phase diagrams in Figure 2. Physical solutions J ν (x, γ; p j ) correspond to the (hatched) finite intervals labelled I 1 . The (dashed) infinite intervals labelled I 2 are associated to unbounded intensities. In this case, subject to J 0 > 0 in phase diagrams (b), (c), (d), (h), (j) or if J 0 ≥ J 1 in phase diagrams (a) and (g), the poles x p (γ, p j ) are determined by withω as the real half-period of ℘ and M integer. Denoting, for a certain M > 0, x p > 0 as the smallest x p (γ, p j ), for all γ in a certain domain, J is physical if h < x p , otherwise (h > x p ) J is unbounded and cannot be used for evaluating the dispersion relation (see Fig. 3). For solutions J ν (x, γ; p j ) to be physical, J 0 and J (h) , with J (h) according to (13), (15), and (16), must be located in intervals I 1 or I 2 (in one of the PDs {R f , J}) and simultaneously in intervals I 1 of PDs {R s,c , J}). These requirements can be summarized as where the condition x p > h refers only to the intervals I 2 . Conditions (A2)-(A7) are presented in Appendix A, where the phase diagrams of Figure 2 are expressed algebraically in terms of R ν and ∂ J R ν (J) using the roots of R ν (J) = 0. We denote condition (18), subject to (15) and (16), as phase diagram condition (PDC) in the following. In parameter space, the PDC defines regions (referred to as PDC regions) that represent parameter sets corresponding to bounded nonnegative solutions J ν± (x, γ; p j ). With respect to the DR (14), we first note, if the PDC is satisfied, that J f± (h, γ; p j ) is continuously differentiable w.r.t. h, γ, p j (due to the PDC, denominators in Eqs. (8) and (9) do not vanish, and, since ℘(h, g 2 , g 3 ) is holomorphic in g 2 (γ, p j ), g 3 (γ, p j ) (see [15], 18.5.1-3), J f± (x, γ; p j ) is continuously differentiable). By varying {h, γ}, p j fixed, and γ, p j satisfying the PDC, both sides of the DR (14) are varying continuously, and a contourplot w.r.t. to {h, γ} can serve as a test for existence of physical solutions {h, γ}. Due to analytic properties mentioned, this procedure also works if h = h 0 is fixed, and a certain parameter p i is variable. This leads to a representation {p i , γ} that describes the dependence of γ on p i . Secondly, conditions (15) and (16) are sufficient for existence of positive roots J (h) . It is necessary to check D ≥ 0 (e.g., by means of a regions plot, p j fixed, in order to select J 0 appropriately for subsequent evaluation). Thirdly, we emphasize that DR (14) must be considered for J f+ and for J f− separately (subject to PDC), leading to two dispersion relations. For simplicity, we use J f+ (x, γ; p j ) in (14), so that (due to Eq. (9)) ∂ x J f+ (x, γ; p j )| x=0 = − R f (J 0 ). Hence, the transmission condition at must be used to select J s+ or J s− appropriately (depending on the signs of µ f and µ s ). As outlined above, the transmission condition at x = h is represented by the DR (14). In general, since ∂xJc±(x,γ;pj )| x=h µc = ± √ Rs(J f+ (h,γ;pj )) µc , the RHS can be positive or negative, depending on the choice of J c− or J c+ , respectively. If a c = 0, the solution J c− (x, γ; p j ) is not consistent with condition (5). Thus, only J c+ (x, γ; p j ) (that remains bounded in this case) can be used, leading to constraint (in (14)) sign Additionally to the PDC, (19) must be satisfied.
Returning to the parameter dependence of solutions γ = γ (h, p j ) of the DR (14), as it stands, it relates γ and h if p j are fixed. In principle, the thickness h does not play a special role. The Weierstrass' function ℘(x; g 2 , g 3 ) in J f± (x, γ; p j ) and J f± (x, γ; p j ) depends analytically on h (if h = 0 modulo real period of ℘(x; g 2 , g 3 )) and g 2 and g 3 (see [15], 18.5.1-18.5.4). Hence ℘(x; g 2 , g 3 ) depends analytically on γ, h, ν , µ ν , a ν , J 0 , and, subject to the PDC, J f± (x, γ; p j ) and J f± (x, γ; p j ) have the same property, so that, in place of h in (14), any of the parameters p i ∈ {γ, h, ν , µ ν , a ν , J 0 } can be considered. The resulting dispersion relations DR (γ, p i ; p j ) depend analytically on γ and on the problem's parameters, and thus can be used to investigate the dependence of γ on any of the parameters p i . Clearly, varying parameters imply varying PDC regions, so that dispersion curves may cross boundaries between different PDC regions. Due to this possibility it is necessary to take all cases, indicated in Figure 2, into account when investigating the dependence of the propagation constant γ (p i ; p j ) on parameters p j . The local behaviour of γ (p i ; p j ) can be used to study the stability of γ w.r.t. p j . A variation p j + δp j of a physical solution γ (p i ; p j ) may lead toγ (p i ; p j + δp j ) that is not consistent with the PDC. It seems suitable to denote such γ (p i ; p j + δp j ) as "basically unstable" with respect to the parameter p j . The other possibilities refer to dispersion curves within the PDC region, where bifurcation points (∂ γ p i = 0, ∂ 2 γ p i = 0) occur. At a bifurcation point γ B the behaviour of the propagation constant γ is described by δγ →0 as δp j → 0 (it depends on the sign of δp j whether two values of γ exist in the vicinity of γ B or one γ on a different branch of the dispersion curve). For the third possibility (δγ → 0 as δp j → 0) it can be shown, due to the continuity of γ = γ (p i ; p j ), that the intensity patterns are continuous. Usually, γ (p i ; p j ) is denoted as "stable" in this case while it is called "bistable" or "unstable" at bifurcation points (see e.g., Sect. 4).
Apart from a positive factor, the time average power flow P tot down the guide is given by For traditional material, γ and P tot are positive. For artificial material, first, γ = 0 is possible (e.g., if s < 0, c < 0), and, second, P is either positive, zero or negative depending on the relative values of three terms. Thus, equation (20) is suitable to study the sign and the strength of the total power flow (see Sect. 4). Following the remarks in the introduction the solution of the problem stated (in Sect. 2) can be summarised as follows.
-If the phase diagram condition (18), subject to (19), are satisfied, solutions γ (p i ; p j ) of the dispersion relation (14) exist. Inserted into equations (6), (8), and (10), they lead to physical intensities J ν± (x, γ; p j ). -Conditions (18) and (19) define certain subspaces in parameter space (see, e.g., A, B in Fig. 3). Branches γ (p i ; p j ) with parameters p j inside or on the boundaries represent parameter dependence of γ (p j ), and hence stability of solutions J ν± (x, γ; p j ) represented by equations (6), (8), and (10). -The analytical representation of the total power flow P according to equation (20), combined with the corresponding parameter dependence of the propagation constant γ can be used to optimize P and to study special effects, e.g., power flow reversal (see Fig. 8).

Numerical example
Before we exploit the foregoing formulae numerically it is necessary to consider the applicability of the model assumptions equations (1) and (2). For traditional material, with positive permittivity and permeability µ, equations (1) and (2) are realistic so that the results of the previous section can be applied [13]. For metamaterial with both negative and µ, this is not self-evident.
First, loss is an intrinsic feature of metamaterials. But with complex and µ, the present approach does not work, so that different methods must be used (see, e.g., [16]).
Second, equations (1) and (2) do not model spatial dispersion at all. Solution of Maxwell's equations with nonlocal (even real) and µ is a non-trivial extension of the foregoing approach that is not in the scope of the present article. In an interesting article Gorlach et al. [17] have presented a theoretical approach for calculating nonlinear susceptibilities of nonlinear discrete metamaterial with both frequency and spatial dispersion. For the structures with negligible spatial dispersion effects expressions for local nonlinear permittivities and susceptibilities are defined (see Eqs. (26), (30), and (31) in [17]) and evaluated numerically. As far as we see, some ad hoc assumptions in equations (17)- (19) [17] lead, admittedly by a very formal argument, to a nonlinear polarization P (ω) = (χ (1) + 3χ (3) E(ω))E(ω), consistent with model equation (1). Further studies of the problem of homogenization of a periodic composite are presented in the literature [18][19][20][21].
Third, it is not clear, whether, in presence of electric and magnetic fields, a magnetic nonlinearity must be taken into account [22]. An extension of the present approach in this direction is (as far as we know) an unsolved problem.  Despite these limitations and concerns there is a widespread use of the above model for metamaterials in the literature [4][5][6][7][8][9][10][11][12].
Finally, the question has to be addressed whether the model assumptions (1) and (2), in particular with f , µ f < 0, are consistent with the Kramers-Kronig Relations (KKRs) (see [23]). As a mathematical theorem, the KKRs are valid if the response of the system is linear, causal, and if the corresponding integrals in the KKRs exist. This means that the response function R(t) is not dependent on the input, which, at time t, should not produce an output for times earlier than t (causality). Finally, R(t) must be square integrable (finite input produces only finite output).
On the other side [26,27], several articles indicate that causality criteria must be used with care. For instance [28], in atomically thin crystals it was found that the KKRs do not hold, confirming a proposition put forward decades ago [29]. As is well known, the KKRs for µ(ω) are not compatible with diamagnetism (µ(0) < 1) and µ > 0 for ω > 0 [30], and in [31] it was pointed out that the general µ(ω) may not satisfy the KKRs, motivating alternative relations to link real and imaginary parts of µ(ω).
The foregoing references indicate that it depends on the physical system and/or the model used whether the KKRs can be applied.
To sum up, there are some limitations related to the model equations (1) and (2), nevertheless, to demonstrate how the approach works, we consider a nonlinear metamaterial film sandwiched between traditional linear cladding and substrate (see Fig. 1) with parameters (for a certain ω) Here, f , µ f have been chosen in accordance to [19] (see Fig. 3) and [10] (see Fig. 1) (it seems that the results in [17] (see Fig. 4) are consistent with this choice). By evaluating the dispersion relation (14), the PDC (18), and P according to (20) we obtain results depicted in Figures 3 and 4. All (red) dispersion curves in Figure 3 are consistent with (18) and (19). As outlined in the captions of Figure 3, branches a, b, c, f correspond to physical solutions, branches d, e are associated to intensities J f+ (x, γ; p j ) with a pole inside the film (see Fig. 5b). For point Q 1 , Figures 4 and 5 show the corresponding phase diagrams (reduced by µ 2 ν , consistent with Fig. 2b) and the field pattern {J s− , J f+ , J c+ }, respectively. For the point Q 4 the field pattern is depicted in Figure 6 (x p < h). Returning to point Q 1 , {h = 1, γ = 1.28}, it is stable w.r.t. the thickness h according to the above mentioned definition of stability. Field patterns of points in the close vicinity of Q 1 are barely distinguishable from field pattern in Figure 5. The question is whether this property is conserved if another parameter is variable. Fixing h = 1 and varying, e.g., c in the vicinity of c = 1 (see parameters (21)), the propagation constant γ is varying due to DR (14) as shown in Figure 7. This procedure can be used for any realistic variation of a (realistic) parameter. Next we consider point Q 2 in Figure 3. As end point of branch b it is specific. Solving the DR (14) with γ = √ s µ s = √ c µ c = 1 one gets h Q2 = 1.818. Since, due to the PDCs (18), γ 2 > ν µ ν , ν = s, c, must hold, Q 2 itself does not represent a physical solution of the dispersion relation (J s− , J c+ are singular). Obviously Q 2 is unstable if it is approached from the left on branch b. For h > h Q2 states S on branch c can be stimulated. Approaching h → h Q2 on branch c from the right, additional to states on branch c, states on branch b (in the vicinity of Q 2 ) can be stimulated. Point Q 2 represents a non-physical state of instability (w.r.t. h). The physical branches a, b, c are outside the regions where P < 0 so that total power flow P tot = γP cannot change sign for these solutions of the dispersion relation. Obviously, for branch f this is not the case. Thus it seems interesting to consider point Q 3 (h Q3 = 2.518, γ Q3 = 1.156) where P tot is changing sign if the thickness h is varying in the vicinity of Q 3 . Since J 0 seems to be more susceptible to experimental control than h we evaluate the dispersion relation (14), the phase diagram conditions (18), and (20) w.r.t. J 0 for fixed h = h Q3 . The result is shown in Figure 8. Increasing J 0 up to bifurcation point S Q3 (J 0 = 1, γ = γ Q3 ), the propagation constant γ jumps to state S on branch g, connected with a jump of the total power flow. To sum up, if the basic assumptions outlined in Sections 1 and 2 are realistic for metamaterial, and if, in particular, metamaterial with parameters according to (21) exist, the results of the present section might have a certain physical significance.
The foregoing example indicates that -the parameter dependence and thus stability behaviour (with respect to material parameters and J 0 ) can be investigated on the basis of PDC, (19), and the compact representation (14) of the DR. -the results of Section 3 are suitable to study the dependence of the total power flow P tot = γP on certain parameters. In particular, if µ ν a ν ≥ 0, ν = s,c, due to (18), vanishing propagation constant γ is impossible. As depicted in Figure 3 (point Q 3 ), P can change sign. Hence power flow reversal is possible.

Conclusions and comments
For a planar waveguide structure, with permittivity and permeability µ according to equations (1) and (2), respectively, we presented the intensities J(x, γ; p j ) of travelling waves (3) as particular solutions of Maxwell's equations (see Eqs. (6), (8), and (9)). We derived the dispersion relation by applying the transmission conditions to solutions J(x, γ; p j ) (see Eq. (14)). For physical solutions J(x, γ; p j ), we presented, using a phase diagram approach, conditions (see (18)) that are suitable to ensure the existence of solutions γ = γ (p i ; p j ) of the DR and to describe the parameter dependence of γ = γ (p i ; p j ), J ν (x, γ; p j ), and P tot = γP . In particular, we investigated power flow reversal.
For traditional and µ in equations (1) and (2) the foregoing conclusions represent an approach that works [13]. As outlined in Section 4, for metamaterial the results must be considered with a certain reservation. Waveguiding in metamaterials with a nonlocal response, nonnegliglible dissipation, and higher-order dispersion cannot be described by the present approach. Nevertheless, the growth of the "Metamaterial Tree of Knowledge" [32] cannot be foreseen, so that the model assumptions of the present approach might become less problematic in the future.
As mentioned, parameter regions where solutions γ of the dispersion relation exist, and the dependence of γ on the parameters have not been studied at all in the related literature (to the best of our knowledge). It seems that these issues are interesting in themselves, and, needless to say, are important for practical applications as exemplified in Section 4.