Eigenvalue Bifurcation in Doubly Nonlinear Problems with an Application to Surface Plasmon Polaritons

We consider a class of generally non-self-adjoint eigenvalue problems which are nonlinear in the solution as well as in the eigenvalue parameter ("doubly"nonlinear). We prove a bifurcation result from simple isolated eigenvalues of the linear problem using a Lyapunov-Schmidt reduction and provide an expansion of both the nonlinear eigenvalue and the solution. We further prove that if the linear eigenvalue is real and the nonlinear problem $\mathcal P\mathcal T$-symmetric, then the bifurcating nonlinear eigenvalue remains real. These general results are then applied in the context of surface plasmon polaritons (SPPs), i.e. localized solutions for the nonlinear Maxwell's equations in the presence of one or more interfaces between dielectric and metal layers. We obtain the existence of transverse electric SPPs in certain $\mathcal P\mathcal T$-symmetric configurations.


Introduction
We study the nonlinear problem where A : L 2 (R d , C) ⊃ D(A) → L 2 (R d , C) is a densely defined, closed (possibly non-self-adjoint) operator with a non-empty resolvent set. Throughout the paper the space L 2 (R d , C) and all other function spaces are complex vector spaces, i.e. defined over the complex field C. The potential W is generally nonlinear in the spectral parameter ω and typically complex valued. The function f is nonlinear in both ω and ϕ and is asymptotically equivalent to a monomial near ϕ = 0. Moreover, we suppose that f is Lipschitz continuous in a neighbourhood of an eigen-pair (ω0, ϕ0), where ω0 ∈ C is a simple isolated eigenvalue of L(x, ·). We prove the bifurcation from ω0 using a fixed point argument and a Lyapunov-Schmidt decomposition. Bifurcation from simple eigenvalues is a well studied problem even in the nonselfadjoint case [6,7,15,18]. In particular, bifurcation in complex Banach spaces (as relevant in our problem) is investigated in [7,15] by means of a Lyapunov-Schmidt reduction coupled with topological degree techniques. However, our result includes also an asymptotic expansion of (ω, ϕ) depending on the behaviour of f for small ϕ and for ω near ω0. More precisely, we find a solution (ω, ϕ) ∈ C × D(A) of the form ω = ω0 + εν + ε 1+τ σ, ϕ = ε α ϕ0 + ε α+1 φ + ε α+1+τ ψ, where ε > 0 is small, ω is the spectral parameter, α is related to the degree of homogeneity of f (x, ω, ·) near 0, τ is a positive parameter, and ν, σ ∈ C as well as φ, ψ ∈ D(A) ∩ ϕ * 0 ⊥ are uniquely determined. Moreover, ν is explicit, see (16), and φ satisfies the linear equation (20).
In [11] the bifurcation was proved (and an asymptotic expansion of (ω, ϕ) was provided) for with A as above. This problem clearly has a linear dependence on the spectral parameter ω. The coefficient ε in (2) is the bifurcation parameter and one studies the bifurcation from an eigenvalue ω0 at ε = 0.
As the bifurcation parameter appears explicitly in the equation, the form of the asymptotic expansion of (ω, ϕ) is unique. Note that (2) can be rescaled to (A − ω)ψ = f (ψ) only for the case of homogeneous nonlinearities f . Therefore, our result extends that of [11] to the case of more general nonlinearities f and a nonlinear dependence of both f and W on the spectral parameter. An important application of non-selfadjoint problems which are nonlinear in ω is the propagation of electromagnetic waves in dispersive media, in particular in structures that include a metal. Interfaces of two different media can support localized waves. A typical example is a surface plasmon polariton (SPP) at the interface of a dielectric and a metal, see e.g. [20,19] or, when more layers of dielectrics and/or metals are considered, [26,24,14]. The general case is, of course, described by Maxwell's equations. Assuming the absence of free charges, we have where E and H is the electric and magnetic field respectively, D = D(E ) is the electric displacement field and ǫ0 and µ0 are respectively the permittivity and the permeability of the free space. The displacement field D is generally nonlinear in E and non-local in time. For odd (e.g. Kerr) nonlinearities and a monochromatic field (E , H, D)(x, y, z, t) = (E, H, D)(x, y, z)e iωt + c.c. (with a real frequency ω) a nonlinear eigenvalue problem in (ω, (E, H)) is obtained if higher harmonics are neglected, for details see Sec. 4. Equation (3) as well as the eigenvalue problem have to be accompanied by interface conditions if an interface of two media is present. Assuming that the interface is planar and parallel to the yz-plane, the interface conditions are E2 = E3 = D1 = 0, H = 0, where we define (for the interface located at x = x0), E2 := E2(x → x + 0 ) − E2(x → x − 0 ), etc., see Sec. 4.
A simple example in the cubically nonlinear case is obtained in structures independent of the y, zvariables by choosing the transverse electric (TE) ansatz E(x, y, z) = (0, 0, ϕ(x)) T e iky (5) with k ∈ R. It leads to the scalar nonlinear problem with functions W, Γ : R 2 → C. The interface conditions here boil down to the continuity condition on ϕ and ϕ ′ , see Section 4.
The bifurcation result provides a curve ε → ω(ε), ϕ(ε) with ω(0) = ω0 and ϕ(0) = 0. Even if ω0 is real, the curve can lie in C \ R for all ε = 0. In order for ω to correspond to the (real) frequency of an electromagnetic field, one needs to ensure that the curve lies in R. As we show in Sec. 3, this is possible by restricting the fixed point argument to a symmetric subspace, namely the PT -symmetric subspace. PT -symmetry has been studied extensively in quantum mechanics, see e.g. [4,5]. Recently, a number of physics papers have studied nonlinear PT -symmetric problems from a phenomenological point of view mainly with emphasis on localized solutions, e.g. [21,27,5]. In the context of SPPs, where metals normally lead to a lossy propagation, PT -symmetry has been applied to obtain lossless propagation, see [2,3]. Mathematically, the restriction of a fixed point argument to a PT -symmetric (or more generally antilinearly symmetric) subspace has been used to obtain real nonlinear eigenvalues, see, e.g., [22,11,9]. This article is organized as follows. In Section 2 we state and prove our main bifurcation result (Theorem 2.1). The realness of the nonlinear eigenvalue is ensured in the case of PT -symmetry in Section 3. Applications to SPPs are then given in Section 4, where 2-and 3-layer-configurations are investigated.
Notation: Henceforth, the norm and the inner product in the underlying Hilbert space L 2 (R d , C) will be denoted by · and ·, · respectively. Moreover, · ∞ stands for the L ∞ (R d ) norm and · A for the graph norm, i.e. u A := u + Au . Note that D(A) = D(L) due to assumptions (A2) and (W1).
for some L > 0 and all ω1,2 ∈ B δ (ω0). An example corresponding to the Drude model for metals (see Sec. 4) is with parameters γ, ωp, k ∈ R andχ (3) a bounded function, Lipschitz continuous in ω. This choice will be important for modelling SPPs. Note that, with such potential W , the operator L(·, ω) := ∂ 2 x +W (·, ω) is non-selfadjoint and also nonlinear in ω. Remark 6. Let us discuss the role of the parameter τ in the expansion (8). Clearly, τ determines the accuracy of the expansion given by the first two terms. According to Theorem 2.1, if αβ ≤ 1, then the optimal value is τ = αβ, which is proportional to the difference of the degree of the lowest degree term in f and the next term. Notice also that higher order terms in the nonlinearity do not play any role in the choice of τ . To give an example, consider the nonlinearities in the table below.
As explained in Sec. 4, in the applications of Theorem 2.1 to time harmonic electromagnetic waves, the relevant nonlinearities are odd. In the case of a cubic nonlinearity (α = 1/2) the first correction term in f is of the kind |ϕ| 4 ϕ and therefore we have β = 2 such that we are in the optimal case τ = 1.
The strategy of the proof of Theorem 2.1 may be summarized as follows. We employ a Lyapunov-Schmidt reduction making use of spectral projections (here the simplicity of the eigenvalue ω0 is used) to decompose the problem into a system in which one rather easily determines ν and φ. Then, the rest becomes a system for the unknowns σ and ψ, which will be solved by means of a nested fixed-point argument. In particular, the assumption that ω0 be isolated is exploited to invert the operator L(·, ω0) restricted to ϕ * 0 ⊥ .
In this initial step we reformulate problem (1) with the ansatz in (8) as a system of two equations using the Lyapunov-Schmidt decomposition.
where v collects all other terms of higher-order in ε, namely Comparing now (13) and (14), the terms of order α + 1 in ε match if and only if we take which is well-defined thanks to assumption (Wt). From the rest of (13)- (14) we obtain Let us now apply Q0 to (1). On the one hand we have and on the other hand Therefore, we obtain An expansion of the last term as in (14)- (15) yields Rewriting (16) as and using (18), we obtain Again, imposing that the terms of the lowest-order in ε match, we get a linear equation for φ: Notice that, with our choice of ν, equation (20) is uniquely solvable in Q0D(A) = D(A) ∩ ϕ * 0 ⊥ by the closed range theorem. Indeed, the operator on the left hand side is Fredholm (see [16,Theorem IV.5.28], where the fact that 0 is a simple isolated eigenvalue of L(·, ω0), see (E1)-(E2), is used) and the right hand side is orthogonal to the kernel of the adjoint operator, i.e. to ϕ * 0 , due to (16).
The rest of (19) produces the following equation for (σ, ψ): Fixed Point Argument.
We proceed by a fixed point argument. Note that although a direct fixed point argument for (σ, ψ) is possible, we opt for a nested version, where we first solve for ψ as a function of σ and subsequently solve for σ. This approach is arguably more transparent.
Let us first address equation (21) and write it as a fixed point equation for ψ, exploiting our assumptions on the eigenvalue ω0, which is assumed simple and isolated. This actually means that Q0L(ω0)Q0 : and its norm is bounded by a constant C(ω0): In our nested fixed point argument for (17), (22) we first solve (22) for ψ ∈ Br 2 (0) ⊂ D(A) for all σ ∈ Br 1 (0) ⊂ C fixed with r1 > 0 arbitrary. To this aim we need to show that for each σ ∈ Br 1 (0) there exists r2 > 0 so that 2 Then, having obtained ψ = ψ(σ) ∈ Br 2 (0), we shall solve equation (17) for σ and finally find a suitable r1. Note that the fixed point argument for equation (17) requires the Lipschitz continuity of σ → ψ(σ), which we verify below.
To ensure (i), we need to estimate R(ψ) . The second and the third term in (21) are easy to handle. Henceforth, we track the dependence of all constants on σ and ψ via r1, r2.
using (9). Let us now deal with R4. Inspecting (15), we obtain where h1 is polynomial in r1, linear in ψ and satisfies h1(0, ψ ) = 0. Actually, all terms appearing in (15) are easy to estimate, so here we just briefly justify the one for the integral rest I(x, ω), for later use too. Indeed, using (11) for n = 3, we have for any 0 < r < δ and all ω ∈ B r/2 (ω0), i.e. for all ε > 0 small enough. Recall that M = sup Finally, we need to estimate R1, which involves the nonlinearity. We split it and estimate term by term. First, by (f4) one gets for ε > 0 small enough. In equation (28) we used the embedding for some r * ∈ [0, r2]. For each r2 > 0 there exists ε0 = ε0(r2) such that for all ε ∈ (0, ε0). Next, where we used the fact that the map z → |z| 1 α z is locally Lipschitz as well as estimates of the type Hence from (27) and therefore, from (23)-(25) and (32), that Using (29) and τ ≤ αβ, we may further estimate for ε > 0 small enough. This follows because for ε small enough we have ε 1+τ h1(r1, r2) < C0.

Bifurcation of real nonlinear eigenvalues: the PT -symmetric case
In applications one often starts from a real eigenvalue of the linear problem and seeks a bifurcation branch (ω, ϕ) where the realness of the nonlinear eigenvalue is preserved, i.e. ω ∈ R. This is required also for the applications to SPPs in Sec. 4. However, the solution ω obtained by Theorem 2.1 is a-priori just complex. In this section we provide a symmetric situation in which the realness of ω is preserved in the bifurcation. This is based on [11, Section III], whereby we adapt the analysis for our more general context. For the sake of completeness and since Section 4 is based on these results, we present most of the details here again.
for all x ∈ R d . Moreover, an operator L acting on a Hilbert space with domain D(L) is PT -symmetric when it commutes with B, namely LB = BL in D(L).
The above operator B is in fact the composition of the operator P, the space reflection (parity), and T , the complex conjugation, which corresponds to the time-reversal in quantum mechanics.
Here |E| 2 = E · E andf (ω) is the Fourier-transform of f . Neglecting higher harmonics is a common approach in theoretical studies of weakly nonlinear optical waves [23]. Using (3), we get that both H and D are curl-fields such that the divergence conditions ∇ · D = 0, ∇ · H = 0 hold automatically. Defining c := (ǫ0µ0) −1/2 , in the second order formulation we have ω 2 Recall again only odd nonlinearities are allowed in D when studying time harmonic waves. Even nonlinearities do not produce terms proportional to e iωt .
For the TE-ansatz in (5), with only one nontrivial component, equation (48) reduces to the scalar problem We study layers of x-periodic (including homogeneous) media. When a metal layer is homogeneous, we choose the simplest Drude model of the linear susceptibility in that layer [19] where ωp ∈ R + and γ ∈ R. For dielectric layers we choose a periodic (possibly constant) and generally complexχ (1) (·, ω). The imaginary part ofχ (1) is related to the loss or gain of energy of an electromagnetic wave propagating inside the medium. For most of the materials, and in particular for metals, it is negative, corresponding to a lossy material. However, in active (doped) materials the energy of an electromagnetic wave is amplified (energy gain), and therefore the imaginary part ofχ (1) is positive. Materials with a realχ (1) (·, ω) are called conservative.
To make equation (49) dimensionless, as well as in order to use the physical values of the parameters involved in the numerical study in Section 4.1.2, we introduce the new rescaled spatial variable, frequency, and wave numberx where ωp is the bulk plasma frequency of a prescribed metal layer. Defining thenφ(x) := ϕ(x), we obtain the same equation as (49) but in the tilde variables and with W and Γ respectively replaced bỹ Note that for the Drude model the susceptibility iŝ For the sake of a simpler notation henceforth we will simply write x, ω, γ, k, W , and Γ instead of x,ω,γ,k,W , andΓ.
Due to the presence of material interface(s), solutions of the Maxwell's equations (3) are not smooth. However, they satisfy the interface conditions that the tangential component of E , the normal component of D and the whole vector H be continuous across each interface, see Sec. 33-3 in [13]. For our interfaces parallel to the yz-plane we get (4). For the ansatz (5) the interface conditions reduce to a C 1 -continuity condition on ϕ ϕ = ∂xϕ = 0.
To apply our bifurcation result to (49) with (IFCs), a real, linear eigenvalue is needed. We show that such an eigenvalue exists in some PT -symmetric choices of the layers. SPPs in PT -symmetric structures have been studied in the physics literature before by employing active materials, see e.g. [2,3,14,25]. Nevertheless, we are not aware of a rigorous mathematical existence proof of the asymptotic expansion of nonlinear SPPs in the frequency dependent case.

The linear eigenvalue problem
In order to apply Theorem 2.1, we need to find an eigen-pair (ω0, ϕ0) of the linear problem coupled with IFCs, with a simple and isolated ω0 ∈ C \ {0}. To ensure the realness of the frequency we need, in fact, ω0 ∈ R \ {0} as well as the PT -symmetry of W (·, ω), such that Proposition 3.1 can be applied. We shall see that the existence of a simple and isolated eigenvalue ω0 ∈ R \ {0} strongly depends on the choice of the layers. As we show in Sec. 4.1.1, the choice of two layers (N = 1) of periodic materials with one being conservative and the other a non-conservative homogeneous material (i.e. with a complexχ (1) ) leads to no real eigenvalues ω0. On the other hand, in Sec. 4.1.2 we find two PT -symmetric settings with three homogeneous layers (N = 2) leading to the existence of an isolated simple eigenvalue ω0 ∈ R in (55). These settings are: (active dielectric -conservative Drude metal -lossy dielectric) and a hypothetical setting of (Drude metal with gain -lossless dielectric -lossy Drude metal).

Two periodic layers
We consider first the case of two layers, each being either a periodic metal or a periodic dielectric, where we set the interface at x = 0. Hence where the functions W±(x, ω) = ω 2 (1 +χ (1) ± (x, ω)) − k 2 are periodic in x with periods ν± > 0. The governing linear problem (55) has two linearly independent Bloch wave solutions ψ ± 1,2 on the half line ±x > 0 respectively. The Bloch wave theory for the Hill's equation (55) can be found in [12]. A necessary condition for the existence of an L 2 (R, C)-solution of (55) is Otherwise (if 0 ∈ S), on at least one of the half lines there is no decaying solution. If 0 / ∈ S, the solutions ψ ± 1,2 have the form where Re(λ±) > 0, p ± j (x + ν±) = p ± j (x) for all x ∈ R and both j = 1, 2. If, say, W+ is real, then p + j and λ+ can be chosen real but p + j is generally 2ν+-periodic, see [12]. An L 2 (R)-solution is given by Due to the linearity of (55) the C 1 -matching condition (IFCs) is equivalent to the condition Note that by varying the parameters ω, k ∈ R we get R± = R±(ω, k).
In [10] the case of real, periodic W±, i.e. the one of two periodic conservative materials, was considered and eigenvalues were found by varying k and searching for zeros of R+(k) − R−(k).
In the presence of non-conservative materials the respective potential W is complex. It is easy to see that the single interface of a conservative material (e.g. a classical dielectric) with a real W and non-conservative homogeneous material (e.g. a Drude metal) with a complex W does not support any eigenvalues of (55). Note that this does not contradict the existence of SPPs at such single metal/dielectric interfaces in general because this existence holds for TM-polarisations, see [19]. Without loss of generality we assume that the conservative material is on the half line x > 0, i.e. W+(x, ω) is real. The nonconservative material in x < 0 is homogeneous (described, e.g. by (50)) such that W−(x, ω) = W−(ω) is complex and independent of x. Hence ϕ−(x) = ce λ − x with 3 λ− = W−(ω) ∈ C \ R, such that R− = λ−, and ϕ+(x) = p + 1 (x)e −λ + x with λ+ ∈ R and p + 1 (x) real and 2ν+−periodic, such that Note that p + 1 (0) = 0 implies c = 0 due to (IFCs) such that only the trivial solution ϕ ≡ 0 is produced in that case.
Because R− ∈ C \ R and R+ ∈ R, condition (56) is not satisfied and no eigenvalue exists.

Three homogeneous layers
Next we consider three homogeneous material layers with interfaces at x = 0 and x = d > 0, i.e. a sandwich geometry with two unbounded layers, where η±, η * ∈ C. A localized solution of (55) is possible only if Note that if the semi-infinite layers are conservative, then η+, η− ∈ R and (58) is equivalent to Under assumption (58) we have where µ := √ −W * , λ± := √ −W±, Re(λ±) > 0, and where A, B, C, D ∈ C are constants to be determined. We can normalize such that D = 1. Then the C 1 -matching at x = 0 and x = d is equivalent to This system has the unique solution together with a condition on the parameter d: The term 2πim appears due the fact that z = log(b) + 2πim solves e z = b for any m ∈ Z. Condition (61) means that under assumption (58) there is a width d > 0 supporting a real eigenvalue if and only if dm =dm(ω) is positive for some value of ω ∈ R and some m ∈ Z. If all the layers are homogeneous conservative materials, (61) cannot be satisfied because µ > 0 and the argument of the logarithm on the right hand side of (61) lies in (0, 1) such that Re(dm) < 0 for all m ∈ Z. Next, we consider three sandwich settings, out of which the last two are PT -symmetric. As the numerical evaluation ofdm suggests, both of these apparently lead to the existence of linear eigenvalues ω0 ∈ R and hence to real bifurcating nonlinear eigenvalues ω. The first one of these PT -symmetric cases has been taken from the physics literature [2,3] while the second one corresponds to a hypothetical material.
Case 1: conservative dielectric -Drude metal -conservative dielectric. Note that a Drude metal layer being sandwiched between two conservative dielectrics does not produce a PT -symmetric potential W . Hence, we do not expect a real eigenvalue ω.
We have η± > 0 and η * = − 1 ω 2 +iγω ∈ C. Our numerical study ofdm shows that for any choice of the constants η+, η−, γ ∈ R and m ∈ Z we cannot find a frequency ω for whichdm(ω) ∈ (0, ∞). Figure 1 (a), (b) shows as an example the behaviour of the maps ω → Re(d0(ω)) and ω → Im(d0(ω)) for different choices of γ, η+, η−. Case 2: PT -symmetric (dielectric -metal -dielectric) setting. For the case of a homogeneous Drude metal sandwiched between two homogeneous non-conservative dielectric layers, we have This setting leads to a PT -symmetric W (with respect to Hence, one of the dielectric layers is lossy while the other is active and generates energy gain. Note that in this example the simple transformation ω ′ := ω 2 leads to a linear dependence on the spectral parameter ω ′ . This configuration of a conservative metal sandwiched between a couple of well-prepared active and lossy dielectrics was considered, e.g., in [2,3]. As the active material we consider titanium dioxide (TiO2), with refractive index n = 3.2 + 0.2i, and as the metal we choose silver with the bulk plasma frequency ωp(Ag) = 8, 85 · 10 15 s −1 . We use ωp(Ag) as the rescaling parameter in (51). Recalling the relation W (x) = ω 2 n(x) 2 − k 2 between the refractive index and the potential W , we choose k = 2 and obtain η± = 9.2 ± 1.28i.

Bifurcation of a nonlinear eigenvalue
We aim to apply the bifurcation result of Theorem 2.1 in its symmetric version provided by Proposition 3.1 in both settings given by Case 2 and Case 3 and find a bifurcating branch of solutions to the reduced Maxwell's equation (49). In both cases we choose the cubic susceptibilityχ (3) ≡ 1. We need thus to verify that ω0, given by the numerical discussion in Section 4.1.2 for a fixed suitable layer width d, is a simple isolated eigenvalue in the sense of (E1)-(E2).
Verification of (E1): ω 0 is simple Because in (59)-(60) the constants A, B, C, D are unique (up to normalization of ϕ0), it is clear that ker L(·, ω0) = span ϕ0, so 0 as eigenvalue of L(·, ω0) is geometrically simple. To prove the algebraic simplicity, suppose by contradiction that there exists a Jordan chain associated to ω0. This means, there exists u ∈ D(A) (see (54)) such that Solving (64) explicitly using the variation of constants, one finds (65) In order to belong to D(A), u must satisfy the C 1 -matching at the interfaces x = 0 and x = d. This implies that the constants c1, c2, c3, c4 have to solve the linear system Note that T is singular since by our choice of d in (61). In order to find a contradiction and exclude the existence of a solution u ∈ D(A) of (64), we now prove that b is not orthogonal to the kernel of T T . Standard computations show that ker T T is one-dimensional and given by The scalar product (b, p) then reads After evaluating the integrals for ϕ0 given by (59)-(60), and some algebraic computations in which the identity e 2µd = (µ−λ + )(µ−λ − ) (µ+λ + )(µ+λ − ) is frequently used, we get Next, we check that (b, p) is non-zero for the values of λ±, µ, and d obtained numerically in Section 4.1.2. We obtain (b, p) ≈ −19.38 − 46.36i for Case 2 and (b, p) ≈ 0.82 + 1.58i for Case 3.
By simple computations, one gets g ′ (0) = g(0) µ(0) This contradicts the fact that d > 0 and that λ± are chosen with positive real part. As a consequence, we deduce that ω0 is isolated in the sense of (E2).
Bifurcation diagrams The numerical computations below were obtained using a centered finite difference discretization of fourth order on an equispaced grid and with the condition that ϕ(x) = 0 for x outside the computational interval. The linear eigenfunction ϕ0 was computed using Matlab's built in eigenvalue solver eigs. The computation of the nonlinear solution ϕ was done via the standard Newton's iteration.