Bifurcation and Asymptotics of Cubically Nonlinear Transverse Magnetic Surface Plasmon Polaritons

Linear Maxwell equations for transverse magnetic (TM) polarized fields support single frequency surface plasmon polaritons (SPPs) localized at the interface of a metal and a dielectric. Metals are typically dispersive, i.e. the dielectric function depends on the frequency. We prove the bifurcation of localized SPPs in dispersive media in the presence of a cubic nonlinearity and provide an asymptotic expansion of the solution and the frequency. The problem is reduced to a system of nonlinear differential equations in one spatial dimension by assuming a plane wave dependence in the direction tangential to the (flat) interfaces. The number of interfaces is arbitrary and the nonlinear system is solved in a subspace of functions with the $H^1$-Sobolev regularity in each material layer. The corresponding linear system is an operator pencil in the frequency parameter due to the material dispersion. Because of the TM-polarization the problem cannot be reduced to a scalar equation. The studied bifurcation occurs at a simple isolated eigenvalue of the pencil. For geometries consisting of two or three homogeneous layers we provide explicit conditions on the existence of eigenvalues and on their simpleness and isolatedness. Real frequencies are shown to exist in the nonlinear setting in the case of PT-symmetric materials. We also apply a finite difference numerical method to the nonlinear system and compute bifurcating curves.


Introduction
In this article we study time harmonic electromagnetic waves at one or more interfaces between layers of nonlinear and dispersive media.In applications these are typically layers of dielectric and metallic materials and the waves are referred to as surface plasmon polaritons (SPPs).The underlying model is given by Maxwell's equations with the absence of free charges, i.e.
In fact, for ωI ̸ = 0, the averaged model ( 4) is not defined.In (1), (2) the linear terms are proportional to e ω I t while the nonlinear terms are proportional to e 3ω I t .This is clearly in contrast with the linear problem ( χ(3) = 0), where complex frequencies are allowed.In fact, the classical linear surface plasmons on metallic surfaces have frequencies with Im(ω) < 0 resulting in an exponential decay in time.
In our analysis we will first prove bifurcations of solutions of (6) with complex ω.Second, we show that under the so called PT -symmetry assumption on the functions χ(1,3) (•, ω) bifurcations with real ω can be proved.The PT -symmetry results in a spatial balance of loss and gain of the material.
The concrete form of the functions χ(1,3) will not be important for our results.These material functions are given by measurements for the given material.A classical description of (homogeneous) metals is via the Drude model [3] where ωp ∈ R + and γ ∈ R. In dielectrics, on the other hand, one often uses an ω-independent approximation of χ (1) .Such an approximation is valid in a certain range of operating frequencies.In homogenous dielectrics one then has χ(1) = const.∈ R.
The realness of χ(1) means that the material is conservative, i.e no energy loss or gain is present.Also χ(3) typically takes the form of a rational function, see Chapter 1 in [7], or it is again assumed to be a real constant (for homogenous materials).SPPs have been observed and studied since the 1950s, see e.g [22], [24] and [23].Classical SPPs arise from the interaction between an illuminating wave and the free electrons of a conductor and generate a highly confined electromagnetic field at the interface of a metal and a dielectric.This makes them useful in sensing (see [2], [18] or for example [16]).The strong localization is useful also in applications requiring light propagation in sub-wavelength geometries [4].In addition, the strong field near the interface enhances the nonlinear response of the medium, see [26], [30], and SPPs propagating in the form of solitons or solitary waves have been studied [9], [14].These and similar studies give a numerical or a formal analytical evidence of such solutions of Maxwell's equations.We give an analytical proof of the existence of localized time harmonic SPPs bifurcating for their linear counterparts.In the case of the transverse electric (TE) polarized fields this bifurcation analysis was carried out in [10,11].Here we address the TM-polarization.Note that in the simplest setting of a single interface between a linear homogenous metal and a linear homogenous dielectric, only TM-SPPs exist.
Interfaces of materials can be understood as waveguides.In fact, a classical waveguide is given by a material sandwiched between two other materials.Such a setting with dielectric homogeneous materials was studied in the Maxwell's equations both for the TE and TM polarization with the finite material being Kerr nonlinear in [29], [31] and [32].Nonlinear dispersion relations were derived in these references.In our case the layers need not be homogenous, all layers can be nonlinear and we allow ω-dependence of the material constants.In addition, the number of layers is arbitrary in our bifurcation result.
As explained in Remark 2, in our time-independent nonlinear equation only solutions with real frequencies ω produce solutions of the nonlinear Maxwell's equations.We achieve the realness of ω by working in PT -symmetric materials.PTsymmetry was originally proposed in quantum mechanics in [6].In the context of SPPs it has been used in [1], [5], [10,11], [15], and [17].Mathematically, we reduce the bifurcation problem to one in a PT -symmetric subspace.Such approach was used also in [12] in a general bifurcation problem but with a linear dependence on the bifurcation parameter.In our case the bifurcation parameter is ω and the dependence on it in the SPP case is generally nonlinear.

Main Results
We consider interfaces of two or more media.These are generally inhomogeneous but with a smooth dependence of χ(1,3) (x, ω) on x within each layer.In the case of m material layers we write where x0 = −∞, xm = ∞, and xj−1 < xj for all j = 1, . . ., m.We use the following notation: where the value of n depends on the context.The cases n = 1, 2, 3 appear in the paper.
Let us now fix the functional analytic setting.Working in the Hilbert space L 2 (R, C 3 ), the domain of the operator A is naturally chosen as ) for every ω in the domain of V (x, •).Note that k is fixed throughout the paper and the operator L k is considered as a pencil with respect to the parameter ω.
Our construction of a family of solutions bifurcating from an eigenvalue ω0 uses the Fredholm property of L k (ω0) -in particular the closedness of its range -as well as the algebraic simpleness of ω0.As L k is an operator pencil, some care has to be given to defining the spectrum as well as simpleness and isolatedness of eigenvalues.We proceed analogously to [8], where the second order formulation corresponding to (6) was studied.Note that unlike in [8] we have not included the linear divergence condition ∇ k • (1 + χ(1) (•, ω))u = 0 in the domain D(A).In [8] this was included to obtain ∇ • D = 0 also if ω = 0.In this paper we are interested only in ω near some eigenvalue ω0 ̸ = 0.In fact, including the divergence condition in the definition of D(A) would complicate the nonlinear analysis as the range of L k (ω) would include only divergence free functions.The nonlinearity h(•, ω, u) is, however, generally not divergence (∇ k •) free for u In applications χ1 is often complex valued.As a result L k (ω) is not self-adjoint and the existence of a real linear eigenvalue ω0 (or a real nonlinear eigenvalue ω) cannot be expected.However, we show that for PT symmetric metamaterials (i.e. with a spatial balance of gain and loss) real linear eigenvalues can be obtained.These persist to real nonlinear eigenvalues.
In this way nonlinear transverse magnetic surface plasmons with real frequencies are found.Such surface plasmons are conservative.
The structure of the rest of the paper is as follows.In Section 2 we first define the spectrum for general operator pencils.Next, we derive explicit conditions for the existence of eigenvalues in the cases of two and three homogeneous layers.We also study their simplicity and isolatedness from the rest of the spectrum.A numerical example is provided where eigenvalues are computed and tested for simplicity and isolatedness.Theorems 1.1 and 1.2 are proved in Sections 3 and 4 respectively.In Section 5 we present a finite difference numerical method for the computation of the bifurcating solutions.Numerical results are shown to confirm the asymptotics given by Theorem 1.1 and Theorem 1.2.Finally, the two appendices provide some supporting calculations for the spectral analysis and for the bifurcation proof.

Linear spectral problem
Similarly to [8] we define the spectrum of the operator pencil L k using an additional parameter λ.Specifically, one considers the standard eigenvalue problem Whether ω belongs to the spectrum (or its subset) of L k is defined below by the condition that λ = 0 belongs to the corresponding set for L k (ω) with ω fixed.We first define the resolvent set of the pencil L k by ) is bijective with a bounded inverse} and the spectrum of The point spectrum is defined by The discrete spectrum is defined via i.e. λ = 0 is an isolated eigenvalue of finite algebraic multiplicity of the standard eigenvalue problem (16) (with ω ∈ Ω fixed).
Here note that the algebraic multiplicity of λ as an eigenvalue of L k (ω) is called infinite if its geometric multiplicity, i.e. dim ker(L k (ω)), is infinite or there exists a sequence (un) n∈N 0 of linearly independent elements un ∈ D(A) such that (L k (ω))un+1 = un for all numbers n ∈ N0 with the function u0 ∈ ker(L k (ω)) \ {0}.Otherwise the algebraic multiplicity is called finite.
Finally ω ∈ σp(L k ) is called algebraically simple if λ = 0 is an algebraically simple eigenvalue of L k (ω), i.e. if it is geometrically simple and there is no solution u ∈ D(A) of where v ∈ ker(L k (ω)) \ {0} and such that u and v are linearly independent.
In order for the spectral theory to be meaningful, we need A (and hence also L k ) to be a closed and densely defined operator.We shall prove the closedness and the denseness results for the operator L k (ω).Proof.The denseness of D(A) in L 2 (R, C 3 ) is obvious.We first show that A : D(A) → L 2 (R, C 3 ) is closed.For this we recall that D(A) equipped with the inner product is a Hilbert space, which is shown completely analogously to Lemma 3.1 in [8].We note that the norm generated by ⟨ and A : This implies that (u k ) is a Cauchy sequence in (D(A), ∥ • ∥A) and hence u ∈ D(A).By the continuity of A we then get Au k → Au in L 2 and hence v = Au.
Proof.We have φ0,2, φ0,3 ∈ H 1 (Ij) for each j by the definition of D(A).The first component is given by φ0,1 =

Homogeneous layers
In the case of homogeneous layers, i.e.Vj(•, ω) = Vj(ω) ∈ C, j ∈ {1, . . ., m}, the fundamental system for the linear part of (6) can be found explicitly and the condition for ω to be an eigenvalue reduces effectively to an algebraic equation.
Next, we consider two special cases, namely m = 2 and m = 3.

Two homogenous layers
Without loss of generality we choose the two homogeneous layers where V±(ω) := V (±x > 0, ω) are independent of x.We also define the functions The spectrum of the problem with two layers was analysed in [8] in detail for the second order formulation (the curl curl problem) with the condition ∇ k • (1 + χ(1) (•, ω))u = 0 included in the definition of the domain of A. As explained in Section 1, this definition of D(A) is not suitable for the nonlinear analysis.Our domain D(A) excludes this condition and hence the range of L k is not divergence free.As a result, in the second order formulation the resolvent problem cannot be reduced to a scalar equation since derivatives of L 2 functions would appear, see the proof of Proposition 3.2 in [8].Hence, there is no substantial benefit in using the second order formulation and we stay within the first order formulation.For this we need to redo some of the calculations in [8] with the definition of D(A) as in (10).
For u3 = 0 we note that + (s)e −µ + s ds, lim Using (30), we get with c independent of r, note that Lemma 3.5 of [8] provides an estimate for all the integral terms in ũ± p : + (s)e µ + s ds For In summary, we get ∥ũ∥ . For u1 we get from ( 24) the estimate . This follows, however, from the estimates (33) because we have, for example, The proof of Proposition 2.3 is complete.
The next result determines the set of eigenvalues of L k away from the set Ω0 for any k ∈ R, and shows their simplicity.
As we show next, outside the set Ω0, the eigenvalues of L k are isolated, i.e. the eigenvalue λ = 0 of L k (ω)u = λu is isolated from the rest of the spectrum of this eigenvalue problem.
Remark 5. Together with Proposition 2.4 we conclude that
As explained below, if dm ∈ (0, ∞) for some m ∈ Z and provided and if (42) does not hold, then with d := dm an eigenvalue of the Maxwell operator L k exists.The equation in (41) can then be reformulated as The term 2πim in (41) appears due the fact that z = log(b) + 2πim solves e z = b for any m ∈ Z. Remark 6. Equation ( 43) is equivalent to the dispersion relation (2.28) in [20].
Note that if all the three layers are conservative materials, i.e.V±, V * : Ω → R+, then (41) cannot be satisfied for any real ω because µ±, µ * > 0 and the absolute value of the right hand side in (43) is smaller than one while the left hand side is larger than one for µ * , d > 0.
For a PT -symmetric example setting (metal with gain -dielectric -metal with loss) we compute dm numerically in Example 1 showing that this setting apparently leads to the existence of linear eigenvalues ω0 ∈ R.
Equation ( 59) has a unique solution C− ∈ C if and only if k) .Note that (42) covers indeed all cases in which the factors of the exponential functions in (60) vanish.The cases µ * = −µ+ with c independent of r, we use again (like in the proof of Proposition 2.3) estimates (33) together with the obvious analogy for the bounded interval (0, d), i.e. .
Next, we prove an analogy to Proposition 2.4, i.e., we determine the point spectrum of L k outside Ω0.
The construction of the above solution ψ ∈ D(A) guarantees the geometric simplicity of any eigenvalue in the sense that λ = 0 is a geometrically simple eigenvalue of L k u = λu.
Next we show that if α(ω) ̸ = 0, then the eigenvalue ω ∈ N (k) is also algebraically simple, which by (19) means that the problem (with ψ as in ( 62)) has no solution u ∈ D(A).Let us set r := ψ ∈ L 2 (R, C 3 ).Assuming for a contradiction that u ∈ D(A) is a solution of (65), the variation of constants from the proof of Proposition 2.6 gives the explicit form (48), (49), and (50) with u1 = i V (x,ω) (r1 − iku3).In order for u to belong to D(A), u2,3 must satisfy the continuity at x = 0 and x = d.This implies that the constants C−, C In order to find a contradiction and exclude the existence of a solution u ∈ D(A) of (65), we now prove that b is not orthogonal to the kernel of T T if α(ω) ̸ = 0. Standard computations show that ker T T is one-dimensional and given by ker T T = span p, p := For the scalar product (b, p) C 4 a direct calculation shows that A number of simplifications, using (43) repeatedly, shows that after multiplication with e µ * d ) the right hand side equals α(ω).
Remark 7. One can easily check in (66) that (b, p) C 4 = 0 if ω ∈ O (k) and the first case in (42) is satisfied.We expect (b, p) C 4 to vanish also in the second case of (42).Remark 8.Note that it can be possibly proved that α(ω) ̸ = 0 for all ω ∈ N (k) .We refrain from trying to show this and only check that α(ω) ̸ = 0 numerically for specific material parameters, see Example 1.
Finally, we study isolatedness of eigenvalues in σp(L k ) \ Ω0.Again, we restrict our attention to eigenvalues in N (k) . where Here, we have used the short notation V * ,± := V * ,±(ω) again.In summary (together with Proposition 2.7) we have for all Proof.Analogously to Proposition (2.5) we choose ω ∈ N (k) and need to show the existence of δ > 0 such that Hence, we work with the linear pencil Let us define µ * ,±(λ Based on Proposition 2.6 we have where ) holds and (68) does not hold}, Clearly, because ω / ∈ (M ) for all λ small enough.It remains to be shown that (67) cannot hold if λ is nonzero and small enough.Otherwise there would be a sequence (λj) ⊂ C such that λj → 0 and f (λj) := e 2µ * (λ j )d −g(λj) = 0 ∀j.Due to the C 1 nature of f , this would imply f ′ (0 . A direct calculation shows that iµ * (0)g ′ (0) g(0)(ω+V * (ω)) = β(ω).
Next, we set the width of the middle layer equal to d1(ω) and check the algebraic simplicity and isolatedness of the corresponding eigenvalue ω ∈ (c, 5) (i.e. the x−coordinate of the point on the graph of the real part in Figure 1 (c)).Using Propositions 2.7 and 2.8, we plot in Figure 2 the quantities α(ω) and |d1(ω) − β(ω)| from Propositions 2.7 and 2.8.Clearly, α remains positive on the whole (c, 5) and d1 ̸ = β everywhere except for at most two points, namely ω ≈ 2.15 and ω ≈ 3.3.Hence, the last three eigenvalues at d = 0.7 mentioned above, i.e. {3.09, 3.44, 4.57} are all algebraically simple and isolated.

Bifurcation of nonlinear surface plasmons
The aim of this section is to prove Theorem 1.1.As stated in this theorem, the solution of ( 6) is expanded in where ω0 ∈ C, φ0 ∈ D(A) are fixed by the linear eigenvalue problem (and the normalization ∥φ0∥ = 1) and ν ∈ C, ϕ ∈ D(A) ∩ ⟨φ * 0 ⟩ ⊥ are fixed by ( 14) and (15).Note that ν is well defined because of the assumption (A-T) and that the normalization ⟨φ0, φ * 0 ⟩ = 1 is allowed because ⟨φ0, φ * 0 ⟩ = 0 would imply the solvability of L k (ω0)v = φ0 for v ∈ D(A), which contradicts the algebraic simplicity of ω0.
Applying P0 to equation ( 6), we get On the other hand, Taylor expanding B(ω) in ω0 up to order two (using assumption (A-V)), we have for any r < δ and ω ∈ Br(ω0).Inserting now the expansions of ω and u from (69) produces where with ω = ω0 + εν + ε 2 σ.
The first inner product on the right hand side of (70) is o(ε 2 ) due to the Lipschitz continuity of h, see (78).Comparing now (70) and (71), the terms of order ε 3 2 match if and only if we set which is well-defined thanks to assumption (A-T).From the rest of ( 70)-(71) we obtain Next, we apply Q0 to (6).First, we write The application of Q0 to (6) produces It remains to solve the system (74), (75) for σ and ψ with ϕ ∈ D(A) ∩ ⟨ϕ * 0 ⟩ ⊥ given by (15).
The H 1 regularity on each layer holds for the second and the third component of each element in D(A) by definition of Remark 9.As we are dealing with an ODE problem, the unique solvability in Lemma 3.1 can be shown explicitly using the variation of parameters; in Appendix A we provide the corresponding calculation for the example of two homogenous layers.An analogous calculation is, in principle, possible for N layers.
To prepare for the fixed point argument, let us start with some estimates on R.
Denote by c a positive constant that may vary from line to line below but be independent of ε, r1, and r2 for ε = ε(r1, r2) > 0 small enough.

Next, let us define
where ∥H∥ W 1,∞ (I j ) for a matrix H denotes the maximum of the W 1,∞ (Ij)-norm of all entries in H.We assume w.l.o.g. that M ≥ |ω0|+|δ| in order to estimate the entry iω in B also by M and obtain max j∈{1,...,m} sup ω∈Br Then, for each j ∈ {1, . . ., m} where q3 is a cubic polynomial with no zero degree terms (such that again q3(εr1) ≤ cεr1 for all ε = ε(r1) small enough).
4 Proof of Theorem 1.2 (bifurcation under the PT -symmetry) We prove next Theorem 1.2, i.e. we show that the under the assumption of PT -symmetry of the coefficients and the realness of ω0, the above bifurcation argument can be carried out in the PT -symmetric subspace resulting in real ω and PT -symmetric u.
As we show now, after multiplication with −i the operator L k (•, ω) as well as the nonlinearity h commute with the PT symmetry provided χ(1,3) are PT -symmetric.We have for u = (u1, u2, u3) T , where Lk (x, ω) := −iL k (x, ω) and h(x, ω, u) := −ih(x, ω, u).For reference, we recall Proof.Let u = (u1, u2, u3) T ∈ D(A).Due to the symmetry of V (x, ω) we get Next, we show the PT -symmetry of the nonlinearity h.As χ(3) (x, ω) is PT -symmetric, one, indeed, obtains Lemma 4.1 allows us to restrict the bifurcation problem to the PT -symmetric subspace of D(A) ∩ H1, i.e. to {u ∈ D(A) ∩ H1 : Bu = u} provided the PT -symmetry of χ(1,3) (•, ω) holds for all ω in a real neighbourhood of ω0.The bifurcating eigenvalue is then real and the solution u is PT -symmetric, i.e.Theorem 1.2 holds.This is proved completely analogously to Proposition 3.1 in [10].
Note that the proof uses the PT -symmetry of the eigenfunction φ0.Under the assumption of algebraic simplicity and realness of ω0 the eigenfunction can always be selected PT -symmetric.Indeed, applying B to the eigenvalue equation yields L k (ω0)(Bφ0) = 0 due to the PT -symmetry of L k .As ω0 is simple, we get φ0 = Bφ0 (up to a multiplicative constant).

Numerical Results
We restrict the numerical computations to two layers (m = 2) with the interface at x = 0. We first reduce the nonlinear system (6) with ω ̸ = 0 to a system of two equations for ũ := (u1, u2).For ω ̸ = 0 the third equation in (6) namely produces and the remaining equations become Due to equation (86) the interface conditions (12) become As ( 87)-( 89) is to be solved for (ũ1, ũ2) and ω, an additional constraint is needed.We choose the condition with a given ε > 0. By varying ε in a right neighbourhood of zero we will thus generate a bifurcation curve in accordance with Theorem 1.1.System (87)-( 90) is solved via the Newton iteration in a finite difference discretization.Because the function |ũ| 2 ũ is not complex differentiable, we work in the real variables ũ1,R, ũ2,R, ũ1,I , ũ2,I , ωR, and ωI , where ũj,R = Re(ũj), ũj,I = Im(ũj), j = 1, 2, ωR = Re(ω), and ωI = Im(ω).We rewrite (87)-(90) as a system of six real equations with four real interface conditions, namely as the real and imaginary parts of equations ( 87)-(90).

Numerical Finite Difference Method
The finite difference discretization is implemented on the interval [−L, L] with L > 0 large enough and with homogenous Dirichlet boundary conditions at x = ±L, which are suitable as we search for localized solutions.Next, we explain that while equation ( 88) is to be solved at all grid points excluding x = 0, equation (87) has to be solved also in the limit x → 0+ and x → 0−.The reason is the condition ∇ • D = 0.The divergence condition is satisfied on R \ {0} by solutions of (87), (88) because these equations imply To get ∇ • D = 0 distributionally on R (which is the correct formulation of the divergence condition for u ∈ D(A)), we need to satisfy also D1 = 0, i.e.V ũ1 − ih1 = 0.This follows automatically from the limits x → 0± of (87) and the second interface condition in (89).After discretization it means that we need to solve (87) also in the limits x → 0− and x → 0+.
Choosing ∆x := 2L N +1 with N ∈ 2N+1, we have the N grid points These are altogether 4N + 4 real degrees of freedom.Solving (87) in the finite difference discretization at xj, j ̸ = j * as well as at x1 → 0± and (88) at xj, j ̸ = j * , we get 2N + 2 real equations from (87) and 2N − 2 equations from (88).Condition (90) produces two real equations.After using the interface condition ikũ1 − ũ′ 2 = 0, two degrees of freedom, e.g.U R 1,j * ,+ and U I 1,j * ,+ are eliminated and we get 4N + 2 real degrees of freedom and 4N + 2 real equations.At xj, j ̸ = j * , we use the centered finite difference stencil of second order for ∂x as well as for ∂ 2 x .At x = 0− and x = 0+ we use the one-sided second order stencils.The integral in the normalization condition (90) is approximated using the trapezoidal rule.

Bifurcation Examples
We present here two numerical examples of bifurcation in the case of 2 layers.One of the examples is PT -symmetric and the linear frequency ω0 (hence also the bifurcating frequency ω) is real.The other example is non-symmetric and the frequency is complex.
The nonlinear equations, discretized using the above finite difference scheme, are solved using the Newton iteration.The initial guess for ω = ω0 + εν with 0 < ε ≪ 1 is provided by the asymptotic approximation ε 1/2 φ0, see Theorem 1.1.The bifurcation curve is obtained by a simple parameter continuation in ε.

PT -symmetric example
We choose k = 1 and a case of a PT -symmetric Drude material that is homogenous on each layer x < 0 and x > 0. In detail, with ωp = 0.5 and γ = 0.7.The resulting function V is plotted in Figure 3 (a).We find a real eigenvalue ω0 by determining real elements of N (k) .In the case of (91) condition (22) can, in fact, be solved explicitly.First, we note that (22)  and we get µ1 = ω 2 0 with ω0 ≈ 1.791.The corresponding eigenfunction φ0 is plotted in Figure 4.In the nonlinearity h we choose ϵ0µ 3 0 χ (3) ≡ 1, i.e. χ (3) is real and x− as well as ω−independent.In particular, it is also PT symmetric.
In Figure 3 (b) we plot the bifurcation diagram showing that the bifurcation parameter (i.e.frequency) ω indeed stays real.As expected, the bifurcation curve is tangent to the line given by the asymptotic approximation ω = ω0 + νε with ν from (14).In Figure 3 (c) we show that the convergence of the asymptotic approximation error |ω − ω0 − εν| is indeed approximately quadratic.The solution u (chosen at ω ≈ 1.7167) plotted in 4 satisfies the PT -symmetry.It is clearly close to the linear solution φ0 but not identical.

Non-PT -symmetric example
We present also one example which is not PT -symmetric, namely with ωp = 0.8, γ = 1, and η = 1.Because of the lack of PT -symmetry, the linear eigenvalues are not real.Equation ( 6) makes sense also for non-real ω.However, as explained in Remark 2, the corresponding solution u does not generate a solution of Maxwell's equations.Nevertheless, we compute here the bifurcation from ω0 ∈ C \ R.
Solving (22) numerically, we obtain an eigenvalue ω0 ≈ 0.4679 − i 0.061.In Fig. 5 (a) we plot V given by (92) with ω = ω0.In (b) we show the bifurcation diagram and the first order asymptotic approximation ω = ω0 + νε.The corresponding asymptotic error is plotted in (c) with an observed approximately quadratic convergence, as predicted.The nonlinear solution u at ω ≈ 0.427 − i 0.066 is plotted in (b) together with the linear eigenfunction φ0 normalized to have a similar amplitude to that of u.

Figure 4 :
Figure 4: The nonlinear solution u at ω ≈ 1.7167.Recall that E 1 = u 1 , E 2 = u 2 , H 3 = u 3 .The linear eigenfunction φ 0 (normalized to have a similar amplitude to that of u) is plotted for comparison.