# Bifurcation and Asymptotics of Cubically Nonlinear Transverse Magnetic
Surface Plasmon Polaritons

Tomáš Dohnal∗ and Runan He  Institut für Mathematik, Martin-Luther-Universität
Halle-Wittenberg,
06120 Halle (Saale), Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(August 28, 2024)

###### Abstract

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. 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 $\mathcal{P}\mathcal{T}$-symmetric materials. We also apply a
finite difference numerical method to the nonlinear system and compute
bifurcating curves.

Keywords: Maxwell equations, surface plasmon, Kerr nonlinearity, bifurcation,
operator pencil, PT-symmetry, asymptotic expansion

##  1 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.

$\partial_{t}\mathcal{D}=\nabla\times\mathcal{H},\quad\mu_{0}\partial_{t}\mathcal{H}=-\nabla\times\mathcal{E},\quad\nabla\cdot\mathcal{D}=\nabla\cdot\mathcal{H}=0,$
(1)

where $\mathcal{E}$ and $\mathcal{H}$ are the electric and magnetic field
respectively, $\mathcal{D}$ is the electric displacement field depending on
$\mathcal{E}$ in a Kerr nonlinear and nonlocal relation

$\begin{split}\mathcal{D}(x,y,z,t)&=\epsilon_{0}\mathcal{E}(x,y,z,t)+\epsilon_{0}\int_{\mathbb{R}}\chi^{(1)}\left(x,y,z,t-s\right)\mathcal{E}(x,y,z,s)\mathrm{d}s\\\
&+\epsilon_{0}\int_{\mathbb{R}}\chi^{(3)}\left(x,y,z,t-s_{1},t-s_{2},t-s_{3}\right)\left(\left(\mathcal{E}(x,y,z,s_{2})\cdot\mathcal{E}(x,y,z,s_{3})\right)\mathcal{E}(x,y,z,s_{1})\right)\mathrm{d}s,\end{split}$
(2)

$\chi^{(1)}:\mathbb{R}^{4}\to\mathbb{R}$,
$\chi^{(3)}:\mathbb{R}^{6}\to\mathbb{R}$, $\chi^{(1)}(\cdot,\tau)=0$ for
$\tau<0$, and $\chi^{(3)}(\cdot,\tau_{1},\tau_{2},\tau_{3})=0$ if $\tau_{1}<0$
or $\tau_{2}<0$ or $\tau_{3}<0$. $\epsilon_{0}$ and $\mu_{0}$ are the
permittivity and the permeability of the free space respectively.

For a monochromatic field

$\left(\mathcal{E},\ \mathcal{H},\ \mathcal{D}\right)(x,y,z,t)=\left(E,\ H,\
D\right)(x,y,z)e^{-\mathrm{i}\omega t}+\left(\overline{E},\ \overline{H},\
\overline{D}\right)(x,y,z)e^{\mathrm{i}\omega t},$

with a real frequency $\omega$, one can obtain a nonlinear eigenvalue problem
in $\left(\omega,\ \left(E,\ H\right)\right)$ by neglecting higher harmonics
(terms proportional to $e^{3\mathrm{i}\omega t}$ and $e^{-3\mathrm{i}\omega
t}$), see [25],

$\nabla\times H=-\mathrm{i}\omega D,\quad\nabla\times
E=\mathrm{i}\omega\mu_{0}H,\quad\nabla\cdot D=\nabla\cdot H=0$ (3a)
$\begin{split}D=\leavevmode\nobreak\
&\epsilon_{0}\left(1+\hat{\chi}^{(1)}(x,y,z,\omega)\right)E\\\
+\leavevmode\nobreak\
&\epsilon_{0}\left((\hat{\chi}^{(3)}(x,y,z,-\omega,\omega,\omega)+\hat{\chi}^{(3)}(x,y,z,\omega,-\omega,\omega))|E|^{2}E+\hat{\chi}^{(3)}(x,y,z,\omega,\omega,-\omega)(E\cdot
E)\overline{E}\right).\end{split}$ (3b)

Here $|E|^{2}=E\cdot\overline{E}$ and $\hat{f}$ is the Fourier transform of
$f$ in time, $\hat{f}(\omega):=\int_{\mathbb{R}}f(t)e^{{\rm i}\omega t}\,{\rm
d}t$. Clearly, if $\omega\neq 0$, then all solutions of the first two
equations in (3a) satisfy the last two equations, i.e. the divergence
conditions.

For notational simplicity we assume

$\hat{\chi}^{(3)}(x,y,z,-\omega,\omega,\omega)=\hat{\chi}^{(3)}(x,y,z,\omega,-\omega,\omega)=\hat{\chi}^{(3)}(x,y,z,\omega,\omega,-\omega)=:\hat{\chi}^{(3)}(x,y,z,\omega).$

Then (3b) becomes

$D=\epsilon_{0}\left(1+\hat{\chi}^{(1)}(x,y,z,\omega)\right)E+\epsilon_{0}\hat{\chi}^{(3)}(x,y,z,\omega)\left(2|E|^{2}E+(E\cdot
E)\overline{E}\right).$

Note that the analysis can be carried out in the same way without this
simplification.

######  Remark 1.

The removal of higher harmonics occurs automatically if one uses a time
averaged model for the nonlinear part of the displacement field, see [21, 27,
28].

We consider one-dimensional structures, where $\hat{\chi}^{(1)}$ and
$\hat{\chi}^{(3)}$ are independent of $y$ and $z$, i.e.

$\hat{\chi}^{(j)}(x,y,z,\omega)=\hat{\chi}^{(j)}(x,\omega),\ j=1,3,$

and study the transverse magnetic fields

$\begin{split}E(x,y,z)&=\left(\mu_{0}u_{1}(x),\mu_{0}u_{2}(x),0\right)^{T}e^{\mathrm{i}ky},\\\
H(x,y,z)&=\left(0,0,u_{3}(x)\right)^{T}e^{\mathrm{i}ky}\end{split}$ (4)

with $k\in\mathbb{R}$ fixed. The factor $\mu_{0}$ has been added for
convenience. This leads to the vectorial nonlinear system

$L_{k}(x,\omega)u:=Au-B(x,\omega)u=h(x,\omega,u),\quad
u:=\left(u_{1},u_{2},u_{3}\right)^{T},\quad x\in\mathbb{R},$ (5)

where

$\begin{split}A:=\begin{pmatrix}0\ &0\ &\mathrm{i}k\\\ 0\ &0\
&-\partial_{x}\\\ -\mathrm{i}k\ &\partial_{x}\
&0\end{pmatrix}=\nabla_{k}\times,\quad\nabla_{k}=\left(\partial_{x},\
\mathrm{i}k,\ 0\right)^{T},\qquad\qquad\\\ B(x,\omega):=\begin{pmatrix}{\rm
i}V(x,\omega)\ &0\ &0\\\ 0\ &{\rm i}V(x,\omega)&0\\\ 0\ &0\ &{\rm
i}\omega\end{pmatrix},\quad
V(x,\omega):=-\omega\epsilon_{0}\mu_{0}\left(1+\hat{\chi}^{(1)}(x,\omega)\right),\\\
h(x,\omega,u):=-{\rm
i}\omega\epsilon_{0}\mu_{0}^{3}\hat{\chi}^{(3)}(x,\omega)\left(2|u_{E}|^{2}u_{E}+(u_{E}\cdot
u_{E})\overline{u_{E}}\right),\qquad\qquad\qquad\quad\end{split}$ (6)

and with $u_{E}:=(u_{1},u_{2},0)^{T}$. In this reduction, solutions $u$ of
system (5) satisfy $\nabla_{k}\cdot(B(\cdot,\omega)u+h(\cdot,\omega,u))=0$,
i.e. $\partial_{x}D_{1}+\partial_{y}D_{2}=0$. Due to the $z$-independence of
$E$ and $H$ and because $H_{1}=H_{2}=0$, we obtain the divergence conditions
$\nabla\cdot D=\nabla\cdot H=0$ as expected. The operator $L_{k}(x,\omega)$
and the matrix $B(x,\omega)$ will often be abbreviated via $L_{k}(\omega)$ and
$B(\omega)$.

######  Remark 2.

Note that problem (5) can be considered also for $\omega\in\mathbb{C}$.
However, the resulting field

$\left(\mathcal{E},\ \mathcal{H},\ \mathcal{D}\right)(x,y,z,t)=\left(E,\ H,\
D\right)(x,y,z)e^{-\mathrm{i}\omega t}+\left(\overline{E},\ \overline{H},\
\overline{D}\right)(x,y,z)e^{\mathrm{i}\overline{\omega}t},$

solves the nonlinear Maxwell equations (3) only if $\omega\in\mathbb{R}$. This
is because if $\operatorname{Im}(\omega)\neq 0$, then the linear terms in (3)
are proportional to $e^{\operatorname{Im}(\omega)t}$ while the nonlinear terms
are proportional to $e^{3\operatorname{Im}(\omega)t}$. This is clearly in
contrast with the linear problem ($\hat{\chi}^{(3)}=0$), where complex
frequencies are allowed. In fact, the classical linear surface plasmons on
metallic surfaces have frequencies with $\operatorname{Im}(\omega)<0$
resulting in an exponential decay in time.

In our analysis we will first prove bifurcations of solutions of (5) with
complex $\omega$. Second, we show that under the so called
$\mathcal{P}\mathcal{T}$-symmetry assumption on the functions
$\hat{\chi}^{(1,3)}(\cdot,\omega)$ bifurcations with real $\omega$ can be
proved. The $\mathcal{P}\mathcal{T}$-symmetry results in a spatial balance of
loss and gain of the material.

The concrete form of the functions $\hat{\chi}^{(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]

$\hat{\chi}^{(1)}(x,\omega)=-\frac{\omega_{p}^{2}}{\omega^{2}+{\rm
i}\gamma\omega},$ (7)

where $\omega_{p}\in\mathbb{R}^{+}$ and $\gamma\in\mathbb{R}$. In dielectrics,
on the other hand, one often uses an $\omega$-independent approximation of
$\hat{\chi}^{(1)}$. Such an approximation is valid in a certain range of
operating frequencies. In homogenous dielectrics one then has
$\hat{\chi}^{(1)}=\text{const.}\in\mathbb{R}.$ The realness of
$\hat{\chi}^{(1)}$ means that the material is conservative, i.e no energy loss
or gain is present. Also $\hat{\chi}^{(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 $\omega$-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 $\omega$ produce solutions of the nonlinear
Maxwell’s equations. We achieve the realness of $\omega$ by working in
$\mathcal{P}\mathcal{T}$-symmetric materials.
$\mathcal{P}\mathcal{T}$-symmetry 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
$\mathcal{P}\mathcal{T}$-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 $\omega$ and
the dependence on it in the SPP case is generally nonlinear.

###  1.1 Main Results

We consider interfaces of two or more media. These are generally inhomogeneous
but with a smooth dependence of $\hat{\chi}^{(1,3)}(x,\omega)$ on $x$ within
each layer. In the case of $m$ material layers we write

$\mathbb{R}=\overline{\cup_{j=1}^{m}I_{j}},\ I_{j}:=(x_{j-1},x_{j}),$

where $x_{0}=-\infty,x_{m}=\infty$, and $x_{j-1}<x_{j}$ for all $j=1,\dots,m$.
We use the following notation:

$\langle\cdot,\cdot\rangle:=\langle\cdot,\cdot\rangle_{L^{2}(\mathbb{R},\mathbb{C}^{n})},\quad\|\cdot\|:=\|\cdot\|_{L^{2}(\mathbb{R})},$
$\mathcal{H}^{1}:=\\{f\in L^{2}(\mathbb{R},\mathbb{C}^{n}):f|_{I_{j}}\in
H^{1}(I_{j},\mathbb{C}^{n})\ \forall j=1,\dots,m\\},\ n\in\mathbb{N},$
$\|f\|_{{\mathcal{H}}^{1}}:=\sum_{j=1}^{m}\|f\|_{H^{1}(I_{j})},$

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}(\mathbb{R},\mathbb{C}^{3})$, the domain of the operator $A$ is
naturally chosen as

$\begin{split}D(A):=\left\\{u\in
L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right):\quad\nabla_{k}\times\ u\in
L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right)\right\\}.\end{split}$ (8)

Assuming that $\hat{\chi}_{1}(\cdot,\omega)\in L^{\infty}(\mathbb{R})$, we
clearly have that $B(\cdot,\omega):D(A)\to L^{2}(\mathbb{R},\mathbb{C}^{3})$
is bounded and $D(L_{k}(\omega))=D(A)$ for every $\omega$ in the domain of
$V(x,\cdot)$. The range of $L_{k}(\omega)$ is
$L^{2}(\mathbb{R},\mathbb{C}^{3})$ for every $\omega$ in the domain of
$V(x,\cdot)$.

We note that

$D(A)=\left\\{u\in L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right):\quad u_{2},\
u_{3}\in H^{1}\left(\mathbb{R},\mathbb{C}\right)\right\\}.$ (9)

For the linear problem, it is useful to study the operator $L_{k}(x,\omega)$
in each material layer separately. Clearly, we have

$\displaystyle D(A)=$ $\displaystyle\left\\{u\in
L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right):\ u_{2},\ u_{3}\in
H^{1}\left(I_{j},\mathbb{C}\right)\ \forall j\in\\{1,\dots,m\\}\text{ and
}\eqref{IFCs}\text{ holds}\right\\}.$ (10) $\llbracket
u_{2}\rrbracket=\llbracket u_{3}\rrbracket=0\ \text{ at }x_{j}\ \forall
j\in\\{1,\dots,m-1\\},$ (11)

where

$\llbracket{f}\rrbracket=0\ \text{ at }x_{j}\text{ means
}\lim_{x\to{x_{j}}^{+}}f\left(x\right)=\lim_{x\to x_{j}^{-}}f(x).$

We assume $V\in L^{\infty}(\mathbb{R},\mathbb{C})$ for all
$\omega\in\Omega\subset\mathbb{C}$ and define

$V_{j}(\cdot,\omega):=V(\cdot,\omega)|_{I_{j}},\quad j\in\\{1,\dots,m\\}$

for each $\omega\in\Omega.$ Defining
$\Omega_{j}:=\\{\omega\in\mathbb{C}:V_{j}(\cdot,\omega)\in
L^{\infty}(\mathbb{R}))\\}$, we have

$\Omega=\cap_{j=1}^{m}\Omega_{j}.$

Note that typically $\Omega\neq\mathbb{C}$ as $\hat{\chi}_{1}$ usually has
poles in the $\omega$-variable, see e.g. the Drude model.

The following assumptions are used in Theorems 1.1 and 1.2.

  1. (A-E)

$\omega_{0}\neq 0$ is an algebraically simple, isolated eigenvalue of the
operator pencil $L_{k}$;

  2. (A-T)

$\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle\neq
0$.

There exists $\delta>0$ such that

  1. (A-V)

$V_{j}(x,\cdot):\mathbb{C}\to\mathbb{C}$ is holomorphic on
$B_{\delta}(\omega_{0})\subset\mathbb{C}$ for each $j\in\\{1,\dots,m\\}$ and
almost every $x\in I_{j}$ and

$V_{j}(\cdot,\omega),\frac{1}{V_{j}(\cdot,\omega)},\partial_{\omega}V_{j}(\cdot,\omega),\partial^{2}_{\omega}V_{j}(\cdot,\omega)\in
W^{1,\infty}(I_{j},\mathbb{C})\quad\forall j\in\\{1,\dots,m\\},\omega\in
B_{\delta}(\omega_{0});$

  2. (A-Na)

$\mathbb{C}\ni\omega\mapsto\hat{\chi}^{(3)}(\cdot,\omega)\in
W^{1,\infty}(I_{j},\mathbb{C})$ is Lipschitz continuous on
$B_{\delta}(\omega_{0})$ for each $j\in\\{1,\dots,m\\}$, i.e. there is
$L_{a}>0$ such that

$\max_{j\in\\{1,\dots,m\\}}\|\hat{\chi}^{(3)}(\cdot,\omega_{1})-\hat{\chi}^{(3)}(\cdot,\omega_{2})\|_{W^{1,\infty}(I_{j})}\leq
L_{a}|\omega_{1}-\omega_{2}|$

for all $\omega_{1},\omega_{2}\in B_{\delta}(\omega_{0})$;

  3. (A-Nb)

$\mathbb{C}\ni\omega\mapsto\hat{\chi}^{(3)}(\cdot,\omega)\in
L^{2}(\mathbb{R})$ is Lipschitz continuous on $B_{\delta}(\omega_{0})$, i.e.
there is $L_{b}>0$ such that

$\|\hat{\chi}^{(3)}(\cdot,\omega_{1})-\hat{\chi}^{(3)}(\cdot,\omega_{2})\|\leq
L_{b}|\omega_{1}-\omega_{2}|$

for all $\omega_{1},\omega_{2}\in B_{\delta}(\omega_{0})$.

######  Theorem 1.1.

Let $k\in\mathbb{R}$. Assume (A-E), (A-V), (A-T), and (A-Na). Let
$\varphi_{0}\in D(A)$ be an eigenfunction corresponding to $\omega_{0}$
normalized to $\|\varphi_{0}\|=1$ and $\varphi_{0}^{*}$ the eigenfunction of
the adjoint $L_{k}^{*}$ with the eigenvalue $\overline{\omega_{0}}$ normalized
to $\langle\varphi_{0},\varphi_{0}^{*}\rangle=1.$

Then there is a unique branch of solutions
$(\omega,u)\in\mathbb{C}\times(D(A)\cap\mathcal{H}^{1})$ of (5) bifurcating
from $(\omega_{0},0)$. There exists $\varepsilon_{0}>0$ s.t. for any
$\varepsilon\in(0,\varepsilon_{0})$ the solution $(\omega,u)$ with $\langle
u,\varphi_{0}^{*}\rangle=\varepsilon^{1/2}$ has the form

$\omega=\omega_{0}+\varepsilon\nu+\varepsilon^{2}\sigma,\qquad\quad
u=\varepsilon^{\frac{1}{2}}\varphi_{0}+\varepsilon^{\frac{3}{2}}\phi+\varepsilon^{\frac{5}{2}}\psi,$
(12)

where

$\nu=-\frac{\langle
h(\cdot,\omega_{0},\varphi_{0}),\varphi_{0}^{*}\rangle}{\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle},$
(13)

$\phi$ is a unique solution in
$D(A)\cap{\mathcal{H}}^{1}\cap\langle\varphi_{0}^{*}\rangle^{\perp}$ of

$L_{k}(x,\omega_{0})\phi=h(x,\omega_{0},\varphi_{0})+\nu\partial_{\omega}B(x,\omega_{0})\varphi_{0},$
(14)

$\sigma\in\mathbb{C}$ and $\psi\in
D(A)\cap{\mathcal{H}}^{1}\cap\langle\varphi_{0}^{*}\rangle^{\perp}.$

######  Theorem 1.2.

In the setting of Theorem 1.1 assume in addition $\omega_{0}\in\mathbb{R}$ and
the $\mathcal{P}\mathcal{T}$-symmetry of the material, i.e.
$V(x,\omega)=\overline{V(-x,\omega)}$ and
$\hat{\chi}^{(3)}(x,\omega)=\overline{\hat{\chi}^{(3)}(-x,\omega)}$ for all
$x\in\mathbb{R}$ and all $\omega\in(\omega_{0}-\delta,\omega_{0}+\delta)$ with
some $\delta>0$. Then the bifurcating solution family with
$\varepsilon\in(0,\varepsilon_{0})$ satisfies $\omega\in\mathbb{R}$ and $u$
can be chosen $\mathcal{P}\mathcal{T}$ symmetric, i.e.
$u(x)=\overline{u(-x)}$.

Our construction of a family of solutions bifurcating from an eigenvalue
$\omega_{0}$ uses the Fredholm property of $L_{k}(\omega_{0})$ \- in
particular the closedness of its range - as well as the algebraic simpleness
of $\omega_{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 (5) was studied. Note that unlike in [8] we have not included
the linear divergence condition
$\nabla_{k}\cdot\left((1+\hat{\chi}^{(1)}(\cdot,\omega))u\right)=0$ in the
domain $D(A)$. In [8] this was included to obtain $\nabla\cdot D=0$ also if
$\omega=0$. In this paper we are interested only in $\omega$ near some
eigenvalue $\omega_{0}\neq 0$. In fact, including the divergence condition in
the definition of $D(A)$ would complicate the nonlinear analysis as the range
of $L_{k}(\omega)$ would include only divergence free functions. The
nonlinearity $h(\cdot,\omega,u)$ is, however, generally not divergence
($\nabla_{k}\cdot$) free for $u$ with
$\nabla_{k}\cdot\left((1+\hat{\chi}^{(1)}(\cdot,\omega))u\right)=0$.

In applications $\hat{\chi}_{1}$ is often complex valued. As a result
$L_{k}(\omega)$ is not self-adjoint and the existence of a real linear
eigenvalue $\omega_{0}$ (or a real nonlinear eigenvalue $\omega$) cannot be
expected. However, we show that for $\mathcal{P}\mathcal{T}$ 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.

##  2 Linear spectral problem

Similarly to [8] we define the spectrum of the operator pencil $L_{k}$ using
an additional parameter $\lambda$. Specifically, one considers the standard
eigenvalue problem

$L_{k}(\omega)u=\lambda u.$ (15)

Whether $\omega$ belongs to the spectrum (or its subset) of $L_{k}$ is defined
below by the condition that $\lambda=0$ belongs to the corresponding set for
$L_{k}(\omega)$ with $\omega$ fixed.

We first define the resolvent set of the pencil $L_{k}$ by

$\rho(L_{k}):=\\{\omega\in\Omega:L_{k}(\omega):D(A)\to
L^{2}(\mathbb{R},\mathbb{C}^{3})\text{ is bijective with a bounded
inverse}\\}$

and the spectrum of $L_{k}$ by

$\sigma(L_{k}):=\Omega\setminus\rho(L_{k}).$ (16)

Note that $\sigma(L_{k})=\\{\omega\in\mathbb{C}:0\in\sigma(L_{k}(\omega))\\}$.

The point spectrum is defined by

$\sigma_{p}(L_{k}):=\\{\omega\in\Omega:\exists u\in
D(A)\setminus\\{0\\}:L_{k}(\omega)u=0\\}.$

Elements of $\sigma_{p}(L_{k})$ are called eigenvalues of $L_{k}$.

The discrete spectrum is defined via

$\omega\in\sigma_{d}(L_{k})\ :\Leftrightarrow\
0\in\sigma_{d}(L_{k}(\omega)),\omega\in\Omega,$ (17)

i.e. $\lambda=0$ is an isolated eigenvalue of finite algebraic multiplicity of
the standard eigenvalue problem (15) (with $\omega\in\Omega$ fixed).

Here note that the algebraic multiplicity of $\lambda$ as an eigenvalue of
$L_{k}(\omega)$ is called infinite if its geometric multiplicity, i.e. ${\rm
dim}\,\ker(L_{k}(\omega))$, is infinite or there exists a sequence
$(u_{n})_{n\in\mathbb{N}_{0}}$ of linearly independent elements $u_{n}\in
D(A)$ such that $(L_{k}(\omega))u_{n+1}=u_{n}$ for all numbers
$n\in\mathbb{N}_{0}$ with the function
$u_{0}\in\ker(L_{k}(\omega))\setminus\\{0\\}$. Otherwise the algebraic
multiplicity is called finite.

Finally $\omega\in\sigma_{p}(L_{k})$ is called algebraically simple if
$\lambda=0$ is an algebraically simple eigenvalue of $L_{k}(\omega)$, i.e. if
it is geometrically simple and there is no solution $u\in D(A)$ of

$L_{k}(\omega)u=v,$ (18)

where $v\in\ker(L_{k}(\omega))\setminus\\{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}(\omega)$.

######  Proposition 2.1.

The operators $A:\ D(A)\longrightarrow
L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right)$ and $L_{k}(\omega):\
D(A)\longrightarrow L^{2}\left(\mathbb{R},\mathbb{C}^{3}\right)$ are closed
and densely defined.

###### Proof.

The denseness of $D(A)$ in $L^{2}(\mathbb{R},\mathbb{C}^{3})$ is obvious. We
first show that $A:D(A)\to L^{2}(\mathbb{R},\mathbb{C}^{3})$ is closed. For
this we recall that $D(A)$ equipped with the inner product

$\langle u,\ v\rangle_{A}:=\langle u,\
v\rangle_{L^{2}(\mathbb{R},\mathbb{C}^{3})}+\langle\nabla_{k}\times u,\
\nabla_{k}\times v\rangle_{L^{2}(\mathbb{R},\mathbb{C}^{3})}$ (19)

is a Hilbert space, which is shown completely analogously to Lemma 3.1 in [8].
We note that the norm generated by $\langle\cdot,\cdot\rangle_{A}$ is the
graph norm

$\|v\|_{A}=\left(\|v\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}^{2}+\|\nabla_{k}\times
v\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}^{2}\right)^{1/2}$

and $A:D(A)\to L^{2}(\mathbb{R},\mathbb{C}^{3})$ is bounded if $D(A)$ is
equipped with $\|\cdot\|_{A}$.

Let now $(u_{k})\subset D(A)$ and $(u_{k},Au_{k})\to(u,v)$ in
$L^{2}(\mathbb{R},\mathbb{C}^{3})^{2}$. This implies that $(u_{k})$ is a
Cauchy sequence in $(D(A),\|\cdot\|_{A})$ and hence $u\in D(A)$. By the
continuity of $A$ we then get $Au_{k}\to Au$ in $L^{2}$ and hence $v=Au$.

Since $V(\cdot,\omega)\in L^{\infty}(\mathbb{R})$, we get that
$B(\cdot,\omega):D(A)\to L^{2}(\mathbb{R},\mathbb{C}^{3})$ is bounded, and
conclude that the operator $L_{k}(\cdot,\omega)$ is a closed and densely
defined operator in $L^{2}(\mathbb{R},\mathbb{C}^{3})$. ∎

######  Lemma 2.2.

Let $\omega_{0}$ be an eigenvalue of $L_{k}$ and assume
$\frac{1}{V_{j}(\cdot,\omega_{0})}\in W^{1,\infty}(I_{j},\mathbb{C})$ for all
$j\in\\{1,\dots,m\\}$. Then corresponding eigenfunctions $\varphi_{0}\in
D(A)\setminus\\{0\\}$ satisfy $\varphi_{0}\in{\mathcal{H}}^{1}.$

###### Proof.

We have $\varphi_{0,2},\varphi_{0,3}\in H^{1}(I_{j})$ for each $j$ by the
definition of $D(A)$. The first component is given by
$\varphi_{0,1}=\frac{k\varphi_{0,3}}{V(\cdot,\omega_{0})}$, which is in
$H^{1}(I_{j})$ for each $j$ as $\frac{1}{V_{j}(\cdot,\omega_{0})}\in
W^{1,\infty}(I_{j})$. ∎

###  2.1 Homogeneous layers

In the case of homogeneous layers, i.e.

$V_{j}(\cdot,\omega)=V_{j}(\omega)\in\mathbb{C},\quad j\in\\{1,\dots,m\\},$

the fundamental system for the linear part of (5) can be found explicitly and
the condition for $\omega$ to be an eigenvalue reduces effectively to an
algebraic equation.

We study eigenvalues in the case of homogeneous layers only outside the set

$\Omega_{0}:=\\{\omega\in\Omega:\omega=0\ \text{ or }\ V_{j}(\omega)=0\text{
for some }j\in\\{1,\dots,m\\}\\}.$

This set consists of eigenvalues of infinite multiplicity. Assuming namely
$V_{j}(\omega)=0$, then
$L_{k}\nabla_{k}\varphi=\nabla_{k}\times(\nabla_{k}\varphi)=0$ for each
$\varphi\in C^{\infty}_{c}(I_{j},\mathbb{R})$ (curl of a gradient vanishes and
$B(x,\omega)(v,w,0)^{T}=0$ for each $x\in I_{j},v,w\in L^{2}(\mathbb{R})$).

In contrast, in [8], where $D(A)$ includes the divergence condition,
$\omega\in\Omega_{0}$ is an eigenvalue of infinite multiplicity only if
$1+\hat{\chi}^{(1)}(\cdot,\omega)=0$ on one of the layers.

Next, we consider two special cases, namely $m=2$ and $m=3$.

####  2.1.1 Two homogenous layers

Without loss of generality we choose the two homogeneous layers
$I_{1}:=\mathbb{R}_{-}:=(-\infty,0),I_{2}:=\mathbb{R}_{+}:=(0,\infty)$ (i.e.
with the interface at $x=0$). Hence

$V(x,\omega)=V_{\pm}(\omega)\quad\text{for }\pm x>0,$

where $V_{\pm}(\omega):=V(\pm x>0,\omega)$ are independent of $x$. We also
define the functions

$W(x,\omega):=-\omega V(x,\omega)\ \text{ and }\ W_{\pm}(\omega):=-\omega
V_{\pm}(\omega).$

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
$\nabla_{k}\cdot\left((1+\hat{\chi}^{(1)}(\cdot,\omega))u\right)=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 (9).

For the resolvent set we have

######  Proposition 2.3.

Let $k\in\mathbb{R}$. Then

$\rho(L_{k})\supset\Omega\setminus(M^{(k)}_{+}\cup M^{(k)}_{-}\cup
N^{(k)}\cup\Omega_{0}),$

where

$M^{(k)}_{\pm}:=\\{\omega\in\Omega\setminus\Omega_{0}:W_{\pm}(\omega)\in[k^{2},\infty)\\},$

and

$N^{(k)}:=\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:W_{+}(\omega),W_{-}(\omega)\notin[k^{2},\infty)\text{
and }\eqref{E:ev.cond-1D}\text{ holds}\\}.$ (20)
$\sqrt{k^{2}-W_{+}(\omega)}V_{-}(\omega)+\sqrt{k^{2}-W_{-}(\omega)}V_{+}(\omega)=0.$
(21)

###### Proof.

Let $r\in L^{2}(\mathbb{R},\mathbb{C}^{3})$ and
$\omega\in\Omega\setminus(M^{(k)}_{+}\cup M^{(k)}_{-}\cup
N^{(k)}\cup\Omega_{0})$. We need to show the existence of a unique $u\in D(A)$
such that

$\left.\begin{array}[]{rcl}{\rm i}ku_{3}-{\rm
i}V_{\pm}(\omega)u_{1}&=&r_{1}\\\ -u_{3}^{\prime}-{\rm
i}V_{\pm}(\omega)u_{2}&=&r_{2}\\\ u_{2}^{\prime}-{\rm i}ku_{1}-{\rm i}\omega
u_{3}&=&r_{3}\end{array}\right\\}{\rm\leavevmode\nobreak\
on\leavevmode\nobreak\ }\mathbb{R}_{\pm},$ (22)

where $V(\omega):=V(x,\omega)$ and such that
$\|u\|_{L^{2}(\mathbb{R})^{3}}\leq c\|r\|_{L^{2}(\mathbb{R})^{3}}$ with $c$
independent of $r$. The first equation in (22) implies

$u_{1}=\frac{{\rm i}}{V_{\pm}(\omega)}\left(r_{1}-{\rm i}ku_{3}\right)\
{\rm\leavevmode\nobreak\ on\leavevmode\nobreak\ }\mathbb{R}_{\pm}.$ (23)

Plugging (23) into the third equatios of (22), one obtains

$\left.\begin{array}[]{rcl}u_{3}^{\prime}+{\rm
i}V_{\pm}(\omega)u_{2}&=&-r_{2}\\\ u_{2}^{\prime}+{\rm
i}\frac{W_{\pm}(\omega)-k^{2}}{V_{\pm}(\omega)}u_{3}&=&r_{3}-\frac{k}{V_{\pm}(\omega)}r_{1}\end{array}\right\\}{\rm\leavevmode\nobreak\
on\leavevmode\nobreak\ }\mathbb{R}_{\pm}.$ (24)

The homogeneous version of system (24) (with the solution vector being
$(u_{2},u_{3})^{T}$) has the fundamental solution matrix

$Y_{\pm}(x):=\begin{pmatrix}\mu_{\pm}e^{\mu_{\pm}x}\
&\mu_{\pm}e^{-\mu_{\pm}x}\\\ -{\rm i}V_{\pm}(\omega)e^{\mu_{\pm}x}\ &{\rm
i}V_{\pm}(\omega)e^{-\mu_{\pm}x}\end{pmatrix},$ (25)

where $\mu_{\pm}:=\sqrt{k^{2}-W_{\pm}(\omega)}$. Analogously to [8] (see Lemma
3.5), we have the following particular solutions in
$L^{2}(\mathbb{R}_{\pm},\mathbb{C}^{2})$

$\displaystyle\tilde{u}_{p}^{+}(x)$
$\displaystyle:=\frac{1}{2}Y_{+}(x)\begin{pmatrix}\int_{x}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s\\\ \int_{0}^{x}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm d}s\end{pmatrix}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}e^{\mu_{+}x}\mu_{+}\int_{x}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s+e^{-\mu_{+}x}\mu_{+}\int_{0}^{x}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm d}s\\\
e^{\mu_{+}x}(-{\rm
i}V_{+}(\omega))\int_{x}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s+e^{-\mu_{+}x}{\rm
i}V_{+}(\omega)\int_{0}^{x}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm
d}s\end{pmatrix},\quad x\in\mathbb{R}_{+},$
$\displaystyle\tilde{u}_{p}^{-}(x)$
$\displaystyle:=\frac{1}{2}Y_{-}(x)\begin{pmatrix}\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s\\\ \int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm d}s\end{pmatrix}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}e^{\mu_{-}x}\mu_{-}\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s+e^{-\mu_{-}x}\mu_{-}\int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\\\ e^{\mu_{-}x}(-{\rm
i}V_{-}(\omega))\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s+e^{-\mu_{-}x}{\rm
i}V_{-}(\omega)\int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\end{pmatrix},\quad x\in\mathbb{R}_{-},$

where

$\rho_{\pm}^{(1)}(x):={\rm
i}\frac{r_{2}(x)}{V_{\pm}(\omega)}+\frac{kr_{1}(x)}{V_{\pm}(\omega)\mu_{\pm}}-\frac{r_{3}(x)}{\mu_{\pm}},\quad\rho_{\pm}^{(2)}(x):={\rm
i}\frac{r_{2}(x)}{V_{\pm}(\omega)}-\frac{kr_{1}(x)}{V_{\pm}(\omega)\mu_{\pm}}+\frac{r_{3}(x)}{\mu_{\pm}}.$
(26)

Because of $\operatorname{Re}(\mu_{\pm})>0$ (recall that $\omega\notin
M_{\pm}^{(k)}$) it holds that $\tilde{u}_{p}^{+}\in
L^{2}(\mathbb{R}_{+},\mathbb{C}^{2})$, $\tilde{u}_{p}^{-}\in
L^{2}(\mathbb{R}_{-},\mathbb{C}^{2})$.

The general solution $\tilde{u}=(u_{2},u_{3})\in
L^{2}(\mathbb{R},\mathbb{C}^{2})$ of (24) is given by

$\tilde{u}(x)=\begin{cases}C_{-}e^{\mu_{-}x}\begin{pmatrix}\mu_{-}\\\ -{\rm
i}V_{-}(\omega)\end{pmatrix}+\tilde{u}_{p}^{-}(x),&x<0,\\\
C_{+}e^{-\mu_{+}x}\begin{pmatrix}\mu_{+}\\\ {\rm
i}V_{+}(\omega)\end{pmatrix}+\tilde{u}_{p}^{+}(x),&x>0\end{cases}$ (27)

with $C_{\pm}\in\mathbb{C}$ arbitrary. For $\omega\notin\Omega_{0}$ we have
$V_{\pm}(\omega)\neq 0$ and (23) yields also $u_{1}\in L^{2}(\mathbb{R}).$

For the interface condition $\llbracket u_{2}\rrbracket=0$ (at $x=0$) first
note that

$\displaystyle\lim_{x\to
0^{+}}u_{2}(x)=C_{+}\mu_{+}+\frac{\mu_{+}}{2}\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s,\quad\lim_{x\to
0^{-}}u_{2}(x)=C_{-}\mu_{-}+\frac{\mu_{-}}{2}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s.$ (28)

Hence $\llbracket u_{2}\rrbracket=0$ if and only if

$C_{+}=\frac{\mu_{-}}{\mu_{+}}C_{-}+\frac{\mu_{-}}{2\mu_{+}}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s-\frac{1}{2}\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm d}s.$
(29)

For $\llbracket u_{3}\rrbracket=0$ we note that

$\begin{split}\displaystyle\lim_{x\to 0^{+}}u_{3}(x)&=C_{+}{\rm
i}V_{+}(\omega)-\frac{{\rm
i}V_{+}(\omega)}{2}\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s,\\\ \lim_{x\to 0^{-}}u_{3}(x)&=-C_{-}{\rm i}V_{-}(\omega)+\frac{{\rm
i}V_{-}(\omega)}{2}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s.\end{split}$ (30)

Using (29), we get

$\left(\mu_{-}V_{+}(\omega)+\mu_{+}V_{-}(\omega)\right)C_{-}=\frac{1}{2}\left(\mu_{+}V_{-}(\omega)-\mu_{-}V_{+}(\omega)\right)\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+\mu_{+}V_{+}(\omega)\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s.$ (31)

Hence, a unique $C_{-}$ exists for any $r\in L^{2}(\mathbb{R},\mathbb{C}^{3})$
if and only if

$\mu_{-}V_{+}(\omega)+\mu_{+}V_{-}(\omega)\neq 0,$

i.e. $\omega\notin N^{(k)}$.

To prove $\|u\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}\leq
c\|r\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}$ with $c$ independent of $r$, note
that Lemma 3.5 of [8] provides an estimate for all the integral terms in
$\tilde{u}_{p}^{\pm}$:

$\displaystyle\left\|e^{\mu_{+}\cdot}\int_{\cdot}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s\right\|_{L^{2}(\mathbb{R}_{+})}$
$\displaystyle\leq\frac{c}{\operatorname{Re}(\mu_{+})}\left\|\rho_{+}^{(1)}\right\|_{L^{2}(\mathbb{R}_{+})},$
(32)
$\displaystyle\left\|e^{-\mu_{+}\cdot}\int_{0}^{\cdot}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm
d}s\right\|_{L^{2}(\mathbb{R}_{+})}$
$\displaystyle\leq\frac{c}{\operatorname{Re}(\mu_{+})}\left\|\rho_{+}^{(2)}\right\|_{L^{2}(\mathbb{R}_{+})},$
$\displaystyle\left\|e^{-\mu_{-}\cdot}\int_{-\infty}^{\cdot}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\right\|_{L^{2}(\mathbb{R}_{-})}$
$\displaystyle\leq\frac{c}{\operatorname{Re}(\mu_{-})}\left\|\rho_{-}^{(2)}\right\|_{L^{2}(\mathbb{R}_{-})},$
$\displaystyle\left\|e^{\mu_{-}\cdot}\int_{\cdot}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s\right\|_{L^{2}(\mathbb{R}_{-})}$
$\displaystyle\leq\frac{c}{\operatorname{Re}(\mu_{-})}\left\|\rho_{-}^{(1)}\right\|_{L^{2}(\mathbb{R}_{-})}.$

For $\omega\in\Omega\setminus(M^{(k)}_{+}\cup M^{(k)}_{-}\cup\Omega_{0})$ we
clearly have

$\|\rho_{\pm}^{(1)}\|_{L^{2}(\mathbb{R}_{\pm})}\leq
c\|r\|_{L^{2}(\mathbb{R}_{\pm},\mathbb{C}^{3})}\ \text{ and }\
\|\rho_{\pm}^{(2)}\|_{L^{2}(\mathbb{R}_{\pm})}\leq
c\|r\|_{L^{2}(\mathbb{R}_{\pm},\mathbb{C}^{3})}.$

In summary, we get $\|\tilde{u}\|_{L^{2}(\mathbb{R},\mathbb{C}^{2})}\leq
c\|r\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}$. For $u_{1}$ we get from (23) the
estimate

$\|u_{1}\|_{L^{2}(\mathbb{R})}\leq
c\left(\|r_{1}\|_{L^{2}(\mathbb{R})}+\|u_{3}\|_{L^{2}(\mathbb{R})}\right).$

Altogether we conclude $\|u\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}\leq
c\|r\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}$.

The only missing property for $u\in D(A)$ is $\tilde{u}\in
H^{1}(\mathbb{R}_{\pm},\mathbb{C}^{2})$. This follows, however, from the
estimates (32) because we have, for example,

$u_{2}^{\prime}(x)=\mu_{-}^{2}\left(C_{-}e^{\mu_{-}x}+\frac{1}{2}\left(e^{\mu_{-}x}\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s-e^{-\mu_{-}x}\int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\right)\right)+\frac{\mu_{-}}{2}(\rho_{-}^{(2)}(x)-\rho_{-}^{(1)}(x))$

for $x\in\mathbb{R}_{-}$.

The proof of Proposition 2.3 is complete. ∎

######  Remark 3.

Note that Proposition 2.3 holds also with $\Omega_{0}$ replaced by the
possibly smaller set
$\tilde{\Omega}_{0}:=\\{\omega\in\Omega:V_{j}(\omega)=0\text{ for some
}j\in\\{1,\dots,m\\}\\}$. In other words, $0\in\rho(L_{k})$ if $0\in\Omega$,
$V_{j}(0)\neq 0\ \forall j$, and $k\neq 0$. The inequality $V_{j}(0)\neq 0$ is
possible if $\hat{\chi}^{(1)}$ has a pole at $\omega=0$. If $k=0$, then
$M_{\pm}^{(k)}=\\{\omega\in\Omega\setminus\tilde{\Omega}_{0}:-\omega
V_{\pm}(\omega)\in[0,\infty)\\}$ and hence $0\in M_{+}^{(k)}\cup M_{-}^{(k)}$
as long as $V_{+}$ or $V_{-}$ has no pole at $0$.

######  Remark 4.

Equation (21) is equivalent to the dispersion relation (2.14) in [20].

The next result determines the set of eigenvalues of $L_{k}$ away from the set
$\Omega_{0}$ for any $k\in\mathbb{R}$, and shows their simplicity.

######  Proposition 2.4.

Let $k\in\mathbb{R}$. Then

$\displaystyle\sigma_{p}(L_{k})\setminus\Omega_{0}$ $\displaystyle=N^{(k)}.$
(33)

All eigenvalues in $\sigma_{p}(L_{k})\setminus\Omega_{0}$ are geometrically
simple. An eigenvalue $\omega\in\sigma_{p}(L_{k})\setminus\Omega_{0}$ is
algebraically simple if and only if

$(2k^{2}-W_{+}(\omega))(2k^{2}-W_{-}(\omega))-2k^{4}-(k^{2}-W_{+}(\omega))V_{-}(\omega)^{2}\neq
0.$ (34)

###### Proof.

From the proof of Prop. 2.3 we get that $L^{2}(\mathbb{R})$-solutions of
$L_{k}(\omega)\psi=0$, i.e. of the homogenous version of (22), have the form

$\psi(x)=\begin{cases}c_{-}e^{\mu_{-}x}\left(\begin{smallmatrix}-{\rm i}k\\\
\mu_{-}\\\ -{\rm i}V_{-}(\omega)\end{smallmatrix}\right)&\quad\text{for}\
x<0,\\\ c_{+}e^{-\mu_{+}x}\left(\begin{smallmatrix}{\rm i}k\\\ \mu_{+}\\\ {\rm
i}V_{+}(\omega)\end{smallmatrix}\right)&\quad\text{for}\ x>0\end{cases}$ (35)

with $c_{\pm}\in\mathbb{C}$. For any $c_{\pm}\in\mathbb{C}$ we clearly have
$\psi|_{\mathbb{R}_{\pm}}\in H^{1}(\mathbb{R}_{\pm},\mathbb{C}^{3})$. The
interface conditions (11) hold if and only if

$c_{-}\mu_{-}=c_{+}\mu_{+}\ \text{ and }\
-c_{-}V_{-}(\omega)=c_{+}V_{+}(\omega),$

which is equivalent to $c_{+}=-\frac{V_{-}(\omega)}{V_{+}(\omega)}c_{-}$ and
$V_{+}(\omega)\mu_{-}+V_{-}(\omega)\mu_{+}=0$, i.e. equation (21).

The above unique form of the eigenfunction $\psi$ shows that each eigenvalue
$\omega$ in $\sigma_{p}({\mathcal{L}}_{k})\setminus\Omega_{0}$ is
geometrically simple in the sense that $\lambda=0$ is a geometrically simple
eigenvalue of

$L_{k}(\omega)u=\lambda u$

Finally, we study the algebraic simpleness, i.e. the non-existence of a
solution $u\in D(A)$ of (18) with $v:=\psi$. That means, we consider

$L_{k}(\omega)u=\psi,\quad u\in D(A).$

Assuming for a contradiction that a solution $u$ exists, we first follow the
lines of the proof of Proposition 2.3 with $r:=\psi\in
L^{2}(\mathbb{R},\mathbb{C}^{3})$. Since (21) holds, we get from (31)

$\displaystyle 0$
$\displaystyle=\frac{1}{2}\left(\mu_{+}V_{-}(\omega)-\mu_{-}V_{+}(\omega)\right)\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+\mu_{+}V_{+}(\omega)\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s,$
$\displaystyle=\mu_{+}V_{-}(\omega)\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+\mu_{+}V_{+}(\omega)\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s=:\alpha,$ (36)

with $\rho_{+}^{(1)}$ and $\rho_{-}^{(2)}$ given in (26).

Next, from (35), setting $c_{-}:=1$ (hence
$c_{+}=-V_{-}(\omega)/V_{+}(\omega)$) and using (21) gives

$\psi_{1}(x)=\begin{cases}-{\rm i}ke^{\mu_{-}x},&x<0,\\\ -{\rm
i}k\frac{V_{-}(\omega)}{V_{+}(\omega)}e^{-\mu_{+}x},&x>0,\end{cases}\quad\psi_{2}(x)=\begin{cases}\mu_{-}e^{\mu_{-}x},&x<0,\\\
\mu_{-}e^{-\mu_{+}x},&x>0,\end{cases}\quad\psi_{3}(x)=\begin{cases}-{\rm
i}V_{-}(\omega)e^{\mu_{-}x},&x<0,\\\ -{\rm
i}V_{-}(\omega)e^{-\mu_{+}x},&x>0\end{cases}$

and we obtain

$\displaystyle\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm d}s=$
$\displaystyle\frac{{\rm
i}}{2\mu_{+}}\left(\frac{\mu_{-}}{V_{+}}-\frac{k^{2}V_{-}}{V_{+}^{2}\mu_{+}}+\frac{V_{-}}{\mu_{+}}\right),$
$\displaystyle\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm d}s=$
$\displaystyle\frac{{\rm
i}}{2\mu_{-}}\left(\frac{\mu_{-}}{V_{-}}+\frac{k^{2}}{V_{-}\mu_{-}}-\frac{V_{-}}{\mu_{-}}\right).$

Using (21), i.e. $\mu_{-}V_{+}+\mu_{+}V_{-}=0$, (36) implies

$\displaystyle-2{\rm i}\frac{\alpha}{V_{+}}=$
$\displaystyle\frac{\mu_{-}}{V_{+}}-\frac{k^{2}V_{-}}{V_{+}^{2}\mu_{+}}+\frac{V_{-}}{\mu_{+}}-\frac{\mu_{-}}{V_{-}}-\frac{k^{2}}{V_{-}\mu_{-}}+\frac{V_{-}}{\mu_{-}}$
$\displaystyle=$
$\displaystyle\mu_{-}\frac{V_{-}-V_{+}}{V_{+}V_{-}}+k^{2}\frac{V_{-}\mu_{-}^{2}-V_{+}\mu_{+}^{2}}{V_{+}\mu_{+}^{2}V_{-}\mu_{-}}+V_{-}\frac{\mu_{+}+\mu_{-}}{\mu_{+}\mu_{-}}$
$\displaystyle=$
$\displaystyle\frac{1}{\mu_{+}^{2}\mu_{-}V_{+}V_{-}}\left[(V_{-}-V_{+})(\mu_{-}^{2}\mu_{+}^{2}+k^{2}(\mu_{+}^{2}+\mu_{-}^{2}))+k^{2}(\mu_{-}^{2}V_{+}-\mu_{+}^{2}V_{-})+\mu_{+}V_{+}V_{-}^{2}(\mu_{+}+\mu_{-})\right].$

Because $\mu_{-}^{2}V_{+}-\mu_{+}^{2}V_{-}=k^{2}(V_{+}-V_{-})$ and
$\mu_{+}V_{+}V_{-}^{2}(\mu_{+}+\mu_{-})=\mu_{+}^{2}V_{-}^{2}(V_{+}-V_{-})$, we
obtain

$\displaystyle-2{\rm i}\frac{\alpha}{V_{+}}=$
$\displaystyle\frac{V_{-}-V_{+}}{\mu_{+}^{2}\mu_{-}V_{+}V_{-}}\left[(\mu_{+}^{2}+k^{2})(\mu_{-}^{2}+k^{2})-2k^{4}-(k^{2}-W_{+})V_{-}^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{V_{-}-V_{+}}{\mu_{+}^{2}\mu_{-}V_{+}V_{-}}\left[(2k^{2}-W_{+})(2k^{2}-W_{-})-2k^{4}-(k^{2}-W_{+})V_{-}^{2}\right].$

Since $V_{+}(\omega)\neq V_{-}(\omega)$ (due to (21)), the statement follows.
∎

As we show next, outside the set $\Omega_{0}$, the eigenvalues of $L_{k}$ are
isolated, i.e. the eigenvalue $\lambda=0$ of $L_{k}(\omega)u=\lambda u$ is
isolated from the rest of the spectrum of this eigenvalue problem.

######  Proposition 2.5.

Let $k\in\mathbb{R}$. Every eigenvalue in
$\sigma_{p}(L_{k})\setminus\Omega_{0}$ is isolated.

###### Proof.

Let $\omega\in\sigma_{p}(L_{k})\setminus\Omega_{0}$. It is to show that
$\lambda\in\rho(A-B(\omega))$ for all $\lambda\in B_{r}(0)\setminus\\{0\\}$
with $r>0$ small enough.

Analogously to Proposition 2.3 we have (replacing $B(\omega)$ by
$B(\omega)+\lambda I$)

$\rho(L_{k})\supset\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:(\omega-{\rm
i}\lambda)(V_{\pm}(\omega)-{\rm i}\lambda)\notin[k^{2},\infty)\text{ and
}\mu_{-}^{(\lambda)}(V_{+}(\omega)-{\rm
i}\lambda)+\mu_{+}^{(\lambda)}(V_{-}(\omega)-{\rm i}\lambda)\neq 0\\},$

where

$\mu_{\pm}^{(\lambda)}:=\sqrt{k^{2}+(\omega-{\rm
i}\lambda)(V_{\pm}(\omega)-{\rm i}\lambda)}.$

Clearly, as $-\omega V_{\pm}(\omega)=W_{\pm}(\omega)\notin[k^{2},\infty)$
(since $\omega\in N^{(k)}$), we get also $-(\omega-{\rm
i}\lambda)(V_{\pm}(\omega)-{\rm i}\lambda)\notin[k^{2},\infty)$ for all
$|\lambda|$ small enough. It remains to be shown that
$\mu_{-}^{(\lambda)}(V_{+}(\omega)-{\rm
i}\lambda)\neq-\mu_{+}^{(\lambda)}(V_{-}(\omega)-{\rm i}\lambda)$ for all
$|\lambda|$ small enough.

Assuming the equality, we get

$(k^{2}+(\omega-{\rm i}\lambda)(V_{-}-{\rm i}\lambda))(V_{+}-{\rm
i}\lambda)^{2}=(k^{2}+(\omega-{\rm i}\lambda)(V_{+}-{\rm
i}\lambda))(V_{-}-{\rm i}\lambda)^{2},$

which can be simplified to

$k^{2}((V_{-}-{\rm i}\lambda)+(V_{+}-{\rm i}\lambda))=-(\omega-{\rm
i}\lambda)(V_{+}-{\rm i}\lambda)(V_{-}-{\rm i}\lambda)$ (37)

after dividing by $V_{+}-V_{-}$, which is non-zero. Equation (37) is a cubic
equation for $\lambda\in\mathbb{C}$. One of the roots is $\lambda=0$ due to
(21), which implies $k^{2}(V_{+}+V_{-})=-\omega V_{+}V_{-}$. We denote the
other two roots by $\lambda_{1},\lambda_{2}$. Choosing
$r<\min\\{|\lambda_{1}|,|\lambda_{2}|\\}$, the equation does not hold for any
$\lambda\in B_{r}(0)\setminus\\{0\\}$. ∎

######  Remark 5.

Together with Proposition 2.4 we conclude that

$\sigma_{d}(L_{k})\setminus\Omega_{0}=\\{\omega\in N^{(k)}:\eqref{E:simple-
cond-2-lay}\text{ holds}\\}.$ (38)

####  2.1.2 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:

$V(x,\omega)=\begin{cases}V_{-}(\omega)&\text{for $x<0$,}\\\
V_{*}(\omega)&\text{for $x\in(0,d)$,}\\\ V_{+}(\omega)&\text{for
$x>d$,}\end{cases}$ (39)

where $V_{-},V_{*},V_{+}:\Omega\subset\mathbb{C}\to\mathbb{C}$. In other
words, we set $I_{1}:=(-\infty,0),I_{2}:=(0,d),I_{3}:=(d,\infty),$
$V_{1}:=V_{-},V_{2}:=V_{*},$ and $V_{3}:=V_{+}$. We also define
$W(x,\omega):=-\omega V(x,\omega)$, $W_{\pm}(\omega):=-\omega
V_{\pm}(\omega)$, and $W_{*}(\omega):=-\omega V_{*}(\omega)$. The following
conditions (an effective dispersion relation) play an important role in the
analysis of Section 2.1.2,

$\exists m\in\mathbb{Z}:\
d_{m}:=\frac{1}{2\mu_{*}}\left[\log\left(\frac{\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{*}V_{-}(\omega)-\mu_{-}V_{*}(\omega)\right)}{\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)}\right)+2\pi{\rm
i}m\right]\in(0,\infty)$ (40)

and

$\begin{cases}&\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)=0,\
\mu_{*}V_{-}(\omega)-\mu_{-}V_{*}(\omega)=0\\\ &\text{or}\\\
&\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)=0,\
\mu_{*}V_{-}(\omega)+\mu_{-}V_{*}(\omega)=0,\end{cases}$ (41)

where

$\mu_{\pm}:=\sqrt{k^{2}-W_{\pm}(\omega)},\quad\text{and}\quad\mu_{*}:=\sqrt{k^{2}-W_{*}(\omega)}.$

As explained below, if $d_{m}\in(0,\infty)$ for some $m\in\mathbb{Z}$ and
provided

$k^{2}-W_{+}(\omega)\notin(-\infty,0],k^{2}-W_{-}(\omega)\notin(-\infty,0],\
\text{and}\ k^{2}-W_{*}(\omega)\neq 0,$

and if (41) does not hold, then with

$d:=d_{m}$

an eigenvalue of the Maxwell operator $L_{k}$ exists. The equation in (40) can
then be reformulated as

$e^{2\mu_{*}d}=\frac{\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{*}V_{-}(\omega)-\mu_{-}V_{*}(\omega)\right)}{\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)}.$
(42)

The term $2\pi{\rm i}m$ in (40) appears due the fact that $z=\log(b)+2\pi{\rm
i}m$ solves $e^{z}=b$ for any $m\in\mathbb{Z}$.

######  Remark 6.

Equation (42) is equivalent to the dispersion relation (2.28) in [20].

Note that if all the three layers are conservative materials, i.e.
$V_{\pm},V_{*}:\Omega\to\mathbb{R}_{+}$, then (40) cannot be satisfied for any
real $\omega$ because $\mu_{\pm},\mu_{*}>0$ and the absolute value of the
right hand side in (42) is smaller than one while the left hand side is larger
than one for $\mu_{*},d>0$.

For a $\mathcal{P}\mathcal{T}$-symmetric example setting (metal with gain -
dielectric - metal with loss) we compute $d_{m}$ numerically in Example 1
showing that this setting apparently leads to the existence of linear
eigenvalues $\omega_{0}\in\mathbb{R}$.

Next, we introduce important sets for the description of the spectrum for any
$k\in\mathbb{R}$:

$\begin{split}M^{(k)}_{-}:=&\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:W_{-}(\omega)\in[k^{2},\infty)\\},\\\
M^{(k)}_{+}:=&\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:W_{+}(\omega)\in[k^{2},\infty)\\},\end{split}$
(43)

and

$N^{(k)}:=\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:W_{+}(\omega),W_{-}(\omega)\notin[k^{2},\infty),W_{*}(\omega)\neq
k^{2},\eqref{cond_d_cD|M|cD}\text{ holds and }\eqref{cond_d_cD|M|cD-2}\text{
does not hold}\\}.$ (44) $O^{(k)}:=\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:W_{+}(\omega),W_{-}(\omega)\notin[k^{2},\infty),W_{*}(\omega)\neq
k^{2},\text{ and }\eqref{cond_d_cD|M|cD-2}\text{ holds}\\}.$ (45)

For the resolvent set we have

######  Proposition 2.6.

Let $k\in\mathbb{R}$.

$\rho(L_{k})\supset\Omega\setminus(M^{(k)}_{+}\cup M^{(k)}_{-}\cup N^{(k)}\cup
O^{(k)}\cup\Omega_{0})=:{\mathcal{M}}^{(k)}.$

###### Proof.

For the whole proof let $\omega\in\Omega\setminus\Omega_{0}.$ As in the proof
of Proposition 2.3, given $r\in L^{2}(\mathbb{R},\mathbb{C}^{3}),$ we firstly
need to show the existence of a unique $u\in D(A)$ such that

$\left.\begin{array}[]{rcl}u_{3}^{\prime}+{\rm i}V(x,\omega)u_{2}&=&-r_{2}\\\
u_{2}^{\prime}+{\rm
i}\frac{W(x,\omega)-k^{2}}{V(x,\omega)}u_{3}&=&r_{3}-\frac{k}{V(x,\omega)}r_{1}\end{array}\right\\}{\rm\leavevmode\nobreak\
on\leavevmode\nobreak\ }(-\infty,0),\ (0,d){\rm\leavevmode\nobreak\
and\leavevmode\nobreak\ }(d,\infty),$ (46)

and $u_{1}:=\frac{{\rm i}}{V(x,\omega)}\left(r_{1}-{\rm i}ku_{3}\right)$.

Analogously to the proof of Proposition 2.3, $L^{2}(\mathbb{R})$-solutions $u$
(disregarding the interface conditions) exist if and only if
$\operatorname{Re}(\mu_{\pm})>0$, i.e.

$k^{2}-W_{\pm}(\omega)\notin(-\infty,0],$

in other words if and only if $\omega\notin M_{+}^{(k)}\cup M_{-}^{(k)}$. The
corresponding $\tilde{u}:=(u_{2},u_{3})^{T}$ has the form

$\tilde{u}(x)=\begin{cases}C_{-}e^{\mu_{-}x}\begin{pmatrix}\mu_{-}\\\ -{\rm
i}V_{-}(\omega)\end{pmatrix}+\tilde{u}_{p}^{-}(x),&x<0,\\\
C_{*}^{(1)}e^{\mu_{*}x}\begin{pmatrix}\mu_{*}\\\ -{\rm
i}V_{*}(\omega)\end{pmatrix}+C_{*}^{(2)}e^{-\mu_{*}x}\begin{pmatrix}\mu_{*}\\\
{\rm i}V_{*}(\omega)\end{pmatrix}+\tilde{u}_{p}^{*}(x),&x\in(0,d),\\\
C_{+}e^{-\mu_{+}x}\begin{pmatrix}\mu_{+}\\\ {\rm
i}V_{+}(\omega)\end{pmatrix}+\tilde{u}_{p}^{+}(x),\quad&x>d,\end{cases}$ (47)

where $C_{\pm},C_{*}^{(1,2)}\in\mathbb{C}$ are arbitrary,

$\displaystyle\tilde{u}_{p}^{-}(x)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}e^{\mu_{-}x}\mu_{-}\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s+e^{-\mu_{-}x}\mu_{-}\int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\\\ -e^{\mu_{-}x}{\rm
i}V_{-}(\omega)\int_{x}^{0}\rho_{-}^{(1)}(s)e^{-\mu_{-}s}\,{\rm
d}s+e^{-\mu_{-}x}{\rm
i}V_{-}(\omega)\int_{-\infty}^{x}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\end{pmatrix},$ (48) $\displaystyle\tilde{u}_{p}^{*}(x)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}e^{\mu_{*}x}\mu_{*}\int_{x}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s+e^{-\mu_{*}x}\mu_{*}\int_{0}^{x}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm d}s\\\
-e^{\mu_{*}x}{\rm
i}V_{*}(\omega)\int_{x}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s+e^{-\mu_{*}x}{\rm
i}V_{*}(\omega)\int_{0}^{x}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s\end{pmatrix},$ $\displaystyle\tilde{u}_{p}^{+}(x)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}e^{\mu_{+}x}\mu_{+}\int_{x}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s+e^{-\mu_{+}x}\mu_{+}\int_{d}^{x}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm d}s\\\
-e^{\mu_{+}x}{\rm
i}V_{+}(\omega)\int_{x}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s+e^{-\mu_{+}x}{\rm
i}V_{+}(\omega)\int_{d}^{x}\rho_{+}^{(2)}(s)e^{\mu_{+}s}\,{\rm
d}s\end{pmatrix},$

and

$\rho_{\pm,*}^{(1)}(x):={\rm
i}\frac{r_{2}(x)}{V_{\pm,*}(\omega)}+\frac{kr_{1}(x)}{V_{\pm,*}(\omega)\mu_{\pm,*}}-\frac{r_{3}(x)}{\mu_{\pm,*}},\quad\rho_{\pm,*}^{(2)}(x):={\rm
i}\frac{r_{2}(x)}{V_{\pm,*}(\omega)}-\frac{kr_{1}(x)}{V_{\pm,*}(\omega)\mu_{\pm,*}}+\frac{r_{3}(x)}{\mu_{\pm,*}}.$
(49)

The regularity $u_{2},u_{3}\in H^{1}(I_{j})$, $j=1,2,3$, follows just like at
the end of the proof of Proposition 2.3 and $u_{1}\in L^{2}(\mathbb{R})$ holds
thanks to $V_{\pm,*}(\omega)\neq 0$ and $r_{1},u_{3}\in L^{2}(\mathbb{R})$.

To show that $u\in D(A)$, it remains to enforce the jump conditions
$\llbracket u_{2}\rrbracket=\llbracket u_{3}\rrbracket=0$ at $x=0,d$. Since

$\displaystyle\lim_{x\to 0^{+}}u_{2}(x)$
$\displaystyle=\left(C_{*}^{(1)}+C_{*}^{(2)}\right)\mu_{*}+\frac{\mu_{*}}{2}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s,$ (50) $\displaystyle\lim_{x\to 0^{-}}u_{2}(x)$
$\displaystyle=C_{-}\mu_{-}+\frac{\mu_{-}}{2}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s,$

and $V_{*}(\omega)\neq 0$, the condition $\llbracket u_{2}\rrbracket=0$ at
$x=0$ is equivalent to

$\left(C_{*}^{(1)}+C_{*}^{(2)}\right)\mu_{*}V_{*}(\omega)=\mu_{-}V_{*}(\omega)C_{-}+\frac{\mu_{-}V_{*}(\omega)}{2}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s-\frac{\mu_{*}V_{*}(\omega)}{2}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s.$ (51)

Similarly, as $\mu_{*}\neq 0$, the interface condition $\llbracket
u_{3}\rrbracket=0$ (at $x=0$), is equivalent to

$\left(C_{*}^{(2)}-C_{*}^{(1)}\right)\mu_{*}V_{*}(\omega)=-\mu_{*}V_{-}(\omega)C_{-}+\frac{\mu_{*}V_{-}(\omega)}{2}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+\frac{\mu_{*}V_{*}(\omega)}{2}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s.$ (52)

Combining (51) and (52), we get

$\begin{split}\mu_{*}V_{*}(\omega)C_{*}^{(1)}&=\frac{\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)}{2}C_{-}+\frac{\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)}{4}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s\\\
&-\frac{\mu_{*}V_{*}(\omega)}{2}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s,\end{split}$ (53)
$\mu_{*}V_{*}(\omega)C_{*}^{(2)}=\frac{\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)}{2}C_{-}+\frac{\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)}{4}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s.$ (54)

Analogously to the above, the interface conditions $\llbracket
u_{2}\rrbracket=0$ and $\llbracket u_{3}\rrbracket=0$ at $x=d$ are

$\begin{split}e^{-\mu_{+}d}\mu_{+}V_{+}(\omega)C_{+}&=\left(C_{*}^{(1)}e^{\mu_{*}d}+C_{*}^{(2)}e^{-\mu_{*}d}\right)\mu_{*}V_{+}(\omega)+\frac{\mu_{*}V_{+}(\omega)}{2}e^{-\mu_{*}d}\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s\\\
&-\frac{\mu_{+}V_{+}(\omega)}{2}e^{\mu_{+}d}\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s.\end{split}$ (55)

and

$\begin{split}&e^{-\mu_{+}d}\mu_{+}V_{+}(\omega)C_{+}=\left(C_{*}^{(2)}e^{-\mu_{*}d}-C_{*}^{(1)}e^{\mu_{*}d}\right)\mu_{+}V_{*}(\omega)+\frac{\mu_{+}V_{*}(\omega)}{2}e^{-\mu_{*}d}\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s\\\
&+\frac{\mu_{+}V_{+}(\omega)}{2}e^{\mu_{+}d}\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s,\end{split}$ (56)

respectively.

Combining (55) and (56), we get

$\begin{split}&e^{\mu_{*}d}\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\mu_{*}V_{*}(\omega)C_{*}^{(1)}+e^{-\mu_{*}d}\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\mu_{*}V_{*}(\omega)C_{*}^{(2)}\\\
&=\frac{\mu_{+}V_{*}(\omega)-\mu_{*}V_{+}(\omega)}{2}\mu_{*}V_{*}(\omega)e^{-\mu_{*}d}\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s+{\mu_{+}V_{+}(\omega)}\mu_{*}V_{*}(\omega)e^{\mu_{+}d}\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s,\end{split}$ (57)

Plugging (53) and (54) into (57), we obtain

$\displaystyle\left[e^{\mu_{*}d}\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)+e^{-\mu_{*}d}\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)\right)\right]C_{-}$
$\displaystyle=(\mu_{+}V_{*}(\omega)-\mu_{*}V_{+}(\omega))\mu_{*}V_{*}(\omega)e^{-\mu_{*}d}\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s+2{\mu_{+}V_{+}(\omega)}\mu_{*}V_{*}(\omega)e^{\mu_{+}d}\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s$
$\displaystyle+(\mu_{+}V_{*}(\omega)+\mu_{*}V_{+}(\omega))\mu_{*}V_{*}(\omega)e^{\mu_{*}d}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s$
$\displaystyle-\frac{1}{2}\left[e^{\mu_{*}d}\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)\right)\right.$
$\displaystyle\qquad\left.+e^{-\mu_{*}d}\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)\right]\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s.$ (58)

Equation (2.1.2) has a unique solution $C_{-}\in\mathbb{C}$ if and only if

$e^{\mu_{*}d}\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)+e^{-\mu_{*}d}\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)\right)\neq
0,$ (59)

which holds because $\omega\notin N^{(k)}\cup O^{(k)}$. Note that (41) covers
indeed all cases in which the factors of the exponential functions in (59)
vanish. The cases
$\mu_{*}=-\mu_{+}\frac{V_{*}}{V_{+}}=\mu_{+}\frac{V_{*}}{V_{+}}$ and
$\mu_{*}=-\mu_{-}\frac{V_{*}}{V_{-}}=\mu_{-}\frac{V_{*}}{V_{-}}$ are both
impossible because $\mu_{\pm},V_{*}(\omega)\neq 0.$

To prove that $\|u\|_{L^{2}(\mathbb{R}^{,}\mathbb{C}^{3})}\leq
c\|r\|_{L^{2}(\mathbb{R},\mathbb{C}^{3})}$ with $c$ independent of $r$, we use
again (like in the proof of Proposition 2.3) estimates (32) together with the
obvious analogy for the bounded interval $(0,d)$, i.e.

$\displaystyle\left\|e^{\mu_{*}\cdot}\int_{\cdot}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm
d}s\right\|_{L^{2}((0,d))}$ $\displaystyle\leq
c\left\|\rho_{*}^{(1)}\right\|_{L^{2}((0,d))},$
$\displaystyle\left\|e^{-\mu_{*}\cdot}\int_{0}^{\cdot}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s\right\|_{L^{2}((0,d))}$ $\displaystyle\leq
c\left\|\rho_{*}^{(2)}\right\|_{L^{2}((0,d))}.$

∎

Next, we prove an analogy to Proposition 2.4, i.e., we determine the point
spectrum of $L_{k}$ outside $\Omega_{0}$.

######  Proposition 2.7.

Let $k\in\mathbb{R}$.

$\sigma_{p}(L_{k})\setminus\Omega_{0}=N^{(k)}\cup O^{(k)}$ (60)

with $N^{(k)}$ given in (44) and $O^{(k)}$ given in (45). All eigenvalues in
$\sigma_{p}(L_{k})\setminus\Omega_{0}$ are geometrically simple. Eigenvalue
$\omega\in N^{(k)}$ is algebraically simple if

$\displaystyle 0\neq\alpha(\omega):=$
$\displaystyle\frac{\mu_{*}^{2}-V_{*}^{2}+k^{2}}{\mu_{*}^{2}V_{*}}\frac{(V_{*}^{2}\mu_{+}\mu_{-}-\mu_{*}^{2}V_{+}V_{-})(\mu_{-}V_{+}+\mu_{+}V_{-})}{\mu_{*}^{2}V_{-}^{2}-\mu_{-}^{2}V_{*}^{2}}$
$\displaystyle+\frac{\mu_{+}^{2}-V_{+}^{2}+k^{2}}{2\mu_{+}^{2}V_{*}}(3\mu_{+}V_{*}-\mu_{*}V_{+})+(\mu_{*}^{2}V_{+}^{2}-\mu_{+}^{2}V_{*}^{2})\left(\frac{1}{\mu_{-}}\frac{\mu_{-}^{2}-V_{-}^{2}+k^{2}}{\mu_{*}^{2}V_{-}^{2}-\mu_{-}^{2}V_{*}^{2}}-d\frac{\mu_{*}^{2}+V_{*}^{2}-k^{2}}{\mu_{*}^{2}V_{*}^{2}}\right),$

where we recall that
$V_{\pm,*}=V_{\pm,*}(\omega),\mu_{\pm,*}=\mu_{\pm,*}(\omega).$

###### Proof.

Let $\omega\in\Omega\setminus\Omega_{0}$ and $k\in\mathbb{R}$. Applying the
arguments of the proof of Proposition 2.6 to the case $r=0$, we get that
$L^{2}(\mathbb{R})$-solutions $\psi$ are given by

$\displaystyle\begin{pmatrix}\psi_{2}\\\ \psi_{3}\end{pmatrix}(x)$
$\displaystyle=\begin{cases}\tilde{A}\left(\begin{smallmatrix}\mu_{+}\\\ {\rm
i}V_{+}(\omega)\end{smallmatrix}\right)e^{-\mu_{+}x},&\quad x>d,\\\
\tilde{B}\left(\begin{smallmatrix}\mu_{*}\\\ {\rm
i}V_{*}(\omega)\end{smallmatrix}\right)e^{-\mu_{*}x}+\tilde{C}\left(\begin{smallmatrix}\mu_{*}\\\
-{\rm i}V_{*}(\omega)\end{smallmatrix}\right)e^{\mu_{*}x}&\quad 0<x<d,\\\
\tilde{D}\left(\begin{smallmatrix}\mu_{-}\\\ -{\rm
i}V_{-}(\omega)\end{smallmatrix}\right)e^{\mu_{-}x},&\quad x<0,\end{cases}$
(61) $\displaystyle\psi_{1}(x)$ $\displaystyle=\frac{k\psi_{3}(x)}{V(\omega)}$

with free constants $\tilde{A}$, $\tilde{B}$, $\tilde{C}$,
$\tilde{D}\in\mathbb{C}$. The $L^{2}$-property holds if and only if
$\omega\notin M_{+}^{(k)}\cup M_{-}^{(k)}$. Once again, $u_{2,3}\in
H^{1}(I_{j}),j=1,2,3$, holds. We can normalize $\psi$ so that $\tilde{D}=1$.
Indeed, if $\tilde{D}=0$, then using the interface conditions, one can show
that $\tilde{A}=\tilde{B}=0$ and hence $\psi\equiv 0.$

Next, we consider the interface conditions (11). The conditions
$\llbracket\psi_{2}\rrbracket=0$ and $\llbracket\psi_{3}\rrbracket=0$ (at
$x=0$ and $x=d$) yield

$\begin{split}\mu_{-}=(\tilde{B}+\tilde{C})\mu_{*},\quad&V_{-}(\omega)=\left(\tilde{C}-\tilde{B}\right)V_{*}(\omega),\quad\tilde{A}e^{-\mu_{+}d}\mu_{+}=(\tilde{B}e^{-\mu_{*}d}+\tilde{C}e^{\mu_{*}d})\mu_{*},\\\
&\tilde{A}e^{-\mu_{+}d}V_{+}(\omega)=\left(\tilde{B}e^{-\mu_{*}d}-\tilde{C}e^{\mu_{*}d}\right)V_{*}(\omega),\quad\end{split}$
(62)

The first equation in (62) yields the necessary condition $\mu_{*}\neq 0$,
i.e. $W_{*}(\omega)\neq k^{2}$. Otherwise $\mu_{-}=0$ and $\psi_{3}\notin
L^{2}(\mathbb{R},\mathbb{C}).$ The first three conditions in (62) now produce
a unique set of coefficients $\tilde{A},\tilde{B},\tilde{C}$, namely

$\displaystyle\tilde{A}$
$\displaystyle=\frac{e^{\mu_{+}d}}{2}\frac{\mu_{*}}{\mu_{+}}\left(e^{\mu_{*}d}\left(\frac{\mu_{-}}{\mu_{*}}+\frac{V_{-}}{V_{*}}\right)+e^{-\mu_{*}d}\left(\frac{\mu_{-}}{\mu_{*}}-\frac{V_{-}}{V_{*}}\right)\right),$
$\displaystyle\tilde{B}$
$\displaystyle=\frac{1}{2}\left(\frac{\mu_{-}}{\mu_{*}}-\frac{V_{-}}{V_{*}}\right),\quad\tilde{C}=\frac{1}{2}\left(\frac{\mu_{-}}{\mu_{*}}+\frac{V_{-}}{V_{*}}\right).$

The last condition in (62) for a nontrivial solution (i.e. for
$\tilde{B},\tilde{C}\neq 0$) is equivalent to

$e^{\mu_{*}d}\left(\mu_{*}V_{+}(\omega)+\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)+\mu_{*}V_{-}(\omega)\right)+e^{-\mu_{*}d}\left(\mu_{*}V_{+}(\omega)-\mu_{+}V_{*}(\omega)\right)\left(\mu_{-}V_{*}(\omega)-\mu_{*}V_{-}(\omega)\right)=0,$
(63)

i.e. (40) or (41), i.e. $\omega\in N^{(k)}\cup O^{(k)}$. An eigenfunction
$\varphi_{0}$ with $\|\varphi_{0}\|=1$ is, of course,
$\varphi_{0}:=\psi/\|\psi\|$.

The construction of the above solution $\psi\in D(A)$ guarantees the geometric
simplicity of any eigenvalue in the sense that $\lambda=0$ is a geometrically
simple eigenvalue of $L_{k}u=\lambda u$.

Next we show that if $\alpha(\omega)\neq 0$, then the eigenvalue $\omega\in
N^{(k)}$ is also algebraically simple, which by (18) means that the problem

$L_{k}(\omega)u=\psi,$ (64)

(with $\psi$ as in (61)) has no solution $u\in D(A)$. Let us set $r:=\psi\in
L^{2}(\mathbb{R},\mathbb{C}^{3})$. Assuming for a contradiction that $u\in
D(A)$ is a solution of (64), the variation of constants from the proof of
Proposition 2.6 gives the explicit form (47), (48), and (49) with
$u_{1}=\frac{{\rm i}}{V(x,\omega)}\left(r_{1}-{\rm i}ku_{3}\right)$.

In order for $u$ to belong to $D(A)$, $u_{2,3}$ must satisfy the continuity at
$x=0$ and $x=d$. This implies that the constants $C_{-}$, $C_{*}^{(1)}$,
$C_{*}^{(2)}$, $C_{+}$ have to solve the linear system

$T\left(C_{-},C_{*}^{(1)},C_{*}^{(2)},C_{+}\right)^{\mathsf{T}}=b,$

where

$T:=\left(\begin{matrix}-\mu_{-}&\mu_{*}&\mu_{*}&0\\\
V_{-}(\omega)&-V_{*}(\omega)&V_{*}(\omega)&0\\\
0&-\mu_{*}e^{\mu_{*}d}&-\mu_{*}e^{-\mu_{*}d}&\mu_{+}e^{-\mu_{+}d}\\\
0&V_{*}(\omega)e^{\mu_{*}d}&-V_{*}(\omega)e^{-\mu_{*}d}&V_{+}(\omega)e^{-\mu_{+}d}\end{matrix}\right)$

and

$b:=\frac{1}{2}\left(\begin{matrix}\mu_{-}\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s-\mu_{*}\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm d}s\\\
V_{-}(\omega)\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+V_{*}(\omega)\int_{0}^{d}\rho_{*}^{(1)}(s)e^{-\mu_{*}s}\,{\rm d}s\\\
e^{-\mu_{*}d}\mu_{*}\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s-e^{\mu_{+}d}\mu_{+}\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s\\\
e^{-\mu_{*}d}V_{*}(\omega)\int_{0}^{d}\rho_{*}^{(2)}(s)e^{\mu_{*}s}\,{\rm
d}s+e^{\mu_{+}d}V_{+}(\omega)\int_{d}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s\end{matrix}\right).$

Note that $T$ is singular since
$T(1,\tilde{C},\tilde{B},\tilde{A})^{\mathsf{T}}=0$, as dictated ba (62). The
functions $\rho^{(1,2)}_{*,\pm}$ in $b$ are as in (49) with $r=\psi.$

In order to find a contradiction and exclude the existence of a solution $u\in
D(A)$ of (64), we now prove that $b$ is not orthogonal to the kernel of
$\overline{T}^{\mathsf{T}}$ if $\alpha(\omega)\neq 0$. Standard computations
show that $\ker\overline{T}^{\mathsf{T}}$ is one-dimensional and given by

$\ker\overline{T}^{\mathsf{T}}=\operatorname{span}\,p,\quad
p:=\begin{pmatrix}e^{\overline{\mu_{*}}d}\overline{V_{-}}(\overline{\mu_{*}}\overline{V_{+}}+\overline{V_{*}}\overline{\mu_{+}})\\\
e^{\overline{\mu_{*}}d}\overline{\mu_{-}}(\overline{\mu_{*}}\overline{V_{+}}+\overline{V_{*}}\overline{\mu_{+}})\\\
\overline{V_{+}}(\overline{\mu_{*}}\overline{V_{-}}-\overline{\mu_{-}}\overline{V_{*}})\\\
-\overline{\mu_{+}}(\overline{\mu_{*}}\overline{V_{-}}-\overline{\mu_{-}}\overline{V_{*}})\\\
\end{pmatrix}.$

For the scalar product $\left(b,p\right)_{\mathbb{C}^{4}}$ a direct
calculation shows that

$\displaystyle-4{\rm i}\left(b,p\right)_{\mathbb{C}^{4}}=$ $\displaystyle
e^{\mu_{*}d}(\mu_{*}V_{+}+\mu_{+}V_{*})\left(\frac{2}{\mu_{-}}(\mu_{-}^{2}-V_{-}^{2}+k^{2})+\frac{d}{\mu_{*}^{2}V_{*}^{2}}(\mu_{-}^{2}V_{*}^{2}-\mu_{*}^{2}V_{-}^{2})(\mu_{*}^{2}+V_{*}^{2}-k^{2})\right)$
(65)
$\displaystyle+\frac{de^{-\mu_{*}d}}{\mu_{*}^{2}V_{*}^{2}}(\mu_{-}V_{*}-\mu_{*}V_{-})^{2}(\mu_{+}V_{*}-\mu_{*}V_{+})(\mu_{*}^{2}+V_{*}^{2}-k^{2})$
$\displaystyle+\frac{\sinh(\mu_{*}d)}{\mu_{*}^{3}V_{*}^{2}}\left((\mu_{-}V_{*}-\mu_{*}V_{-})^{2}(\mu_{*}V_{+}+\mu_{+}V_{*})+(\mu_{*}^{2}V_{-}^{2}-\mu_{-}^{2}V_{*}^{2})(\mu_{*}V_{+}-\mu_{+}V_{*}))\right)(\mu_{*}^{2}-V_{*}^{2}+k^{2})$
$\displaystyle+\frac{2}{\mu_{+}^{2}V_{*}}(\mu_{-}V_{*}\cosh(\mu_{*}d)+\mu_{*}V_{-}\sinh(\mu_{*}d))(\mu_{-}V_{*}-\mu_{*}V_{-})(\mu_{+}^{2}-V_{+}^{2}+k^{2}).$

A number of simplifications, using (42) repeatedly, shows that after
multiplication with
$\frac{e^{\mu_{*}d}}{2}\frac{\mu_{*}V_{+}+\mu_{+}V_{*}}{(\mu_{*}V_{-}-\mu_{-}V_{*})^{2}}$
(which is non-zero since $\omega\notin O^{(k)}$) the right hand side equals
$\alpha(\omega)$. ∎

######  Remark 7.

One can easily check in (65) that $(b,p)_{\mathbb{C}^{4}}=0$ if $\omega\in
O^{(k)}$ and the first case in (41) is satisfied. We expect
$(b,p)_{\mathbb{C}^{4}}$ to vanish also in the second case of (41).

######  Remark 8.

Note that it can be possibly proved that $\alpha(\omega)\neq 0$ for all
$\omega\in N^{(k)}$. We refrain from trying to show this and only check that
$\alpha(\omega)\neq 0$ numerically for specific material parameters, see
Example 1.

Finally, we study isolatedness of eigenvalues in
$\sigma_{p}(L_{k})\setminus\Omega_{0}$. Again, we restrict our attention to
eigenvalues in $N^{(k)}$.

######  Proposition 2.8.

Let $k\in\mathbb{R}$. Eigenvalue $\omega\in N^{(k)}$ is isolated if and only
if

$d\neq\beta(\omega),$

where

$\displaystyle\beta(\omega):=$
$\displaystyle\frac{1}{(\omega+V_{*})(V_{*}^{2}\mu_{-}^{2}-V_{-}^{2}\mu_{*}^{2})(V_{*}^{2}\mu_{+}^{2}-V_{+}^{2}\mu_{*}^{2})}\left\\{(V_{*}(\omega+V_{*})-2\mu_{*}^{2})(\mu_{*}^{2}V_{+}V_{-}-V_{*}^{2}\mu_{+}\mu_{-})(\mu_{+}V_{-}+\mu_{-}V_{+})\right.$
$\displaystyle\left.-\frac{\mu_{*}^{2}V_{*}}{\mu_{+}\mu_{-}}\left[\mu_{-}V_{+}(\omega+V_{+})(\mu_{*}^{2}V_{-}^{2}-V_{*}^{2}\mu_{-}^{2})+\mu_{+}V_{-}(\omega+V_{-})(\mu_{*}^{2}V_{+}^{2}-V_{*}^{2}\mu_{+}^{2})\right]\right.$
$\displaystyle\left.+2\mu_{*}^{2}V_{*}\left[\mu_{*}^{2}(\mu_{+}V_{-}^{2}+\mu_{-}V_{+}^{2})-V_{*}^{2}\mu_{+}\mu_{-}(\mu_{+}+\mu_{-})\right]\right\\}.$

Here, we have used the short notation $V_{*,\pm}:=V_{*,\pm}(\omega)$ again.

In summary (together with Proposition 2.7) we have for all $k\in\mathbb{R}$

$\\{\omega\in N^{(k)}:\alpha(\omega)\neq 0,\beta(\omega)\neq
d\\}\subset\sigma_{d}(L_{k})\setminus\Omega_{0}.$

###### Proof.

Analogously to Proposition (2.5) we choose $\omega\in N^{(k)}$ and need to
show the existence of $\delta>0$ such that

$\lambda\in B_{\delta}(0)\setminus\\{0\\}\subset\mathbb{C}\ \Rightarrow\
\lambda\in\rho(L_{k}(\omega)).$

Hence, we work with the linear pencil
$L_{k}^{(\lambda)}(\omega):=L_{k}(\omega)-\lambda I=A-(B(\omega)+\lambda I)$
and need to show

$\omega\in\rho(L_{k}^{(\lambda)})\quad\forall\lambda\in
B_{\delta}(0)\setminus\\{0\\}\subset\mathbb{C}.$

Let us define

$\mu_{*,\pm}(\lambda):=\sqrt{k^{2}+(\omega-{\rm
i}\lambda)(V_{*,\pm}(\omega)-{\rm i}\lambda)},\quad
V_{*,\pm}(\omega,\lambda):=V_{*,\pm}(\omega)-{\rm i}\lambda.$

Based on Proposition 2.6 we have

$\rho(L_{k}^{(\lambda)})\supset\Omega\setminus(M^{(k,\lambda)}_{+}\cup
M^{(k,\lambda)}_{-}\cup N^{(k,\lambda)}\cup
O^{(k,\lambda)}\cup\Omega_{0}^{(\lambda)}),$

where

$\begin{split}M^{(k,\lambda)}_{\pm}:=&\leavevmode\nobreak\
\\{\omega\in\Omega\setminus\Omega_{0}:-(\omega-{\rm
i}\lambda)V_{\pm}(\omega,\lambda)\in[k^{2},\infty)\\},\\\
N^{(k,\lambda)}:=&\leavevmode\nobreak\
\\{\omega\in\Omega\setminus(\Omega_{0}\cup M^{(k,\lambda)}_{+}\cup
M^{(k,\lambda)}_{-}):-(\omega-{\rm i}\lambda)V_{*}(\omega,\lambda)\neq
k^{2},\eqref{cond_d_cD|M|cD-lam}\text{ holds and
}\eqref{cond_d_cD|M|cD-2-lam}\text{ does not hold}\\},\\\
O^{(k,\lambda)}:=&\leavevmode\nobreak\
\\{\omega\in\Omega\setminus(\Omega_{0}\cup M^{(k,\lambda)}_{+}\cup
M^{(k,\lambda)}_{-}):\eqref{cond_d_cD|M|cD-2-lam}\text{ holds}\\},\\\
\Omega_{0}^{(\lambda)}:=&\leavevmode\nobreak\ \\{\omega\in\Omega:\omega-{\rm
i}\lambda=0\text{ or }V_{-}(\omega,\lambda)=0\ \text{or}\
V_{*}(\omega,\lambda)=0\ \text{or}\ V_{+}(\omega,\lambda)=0\\},\end{split}$
$e^{2\mu_{*}(\lambda)d}=\frac{\left(\mu_{*}(\lambda)V_{+}(\omega,\lambda)-\mu_{+}(\lambda)V_{*}(\omega,\lambda)\right)\left(\mu_{*}(\lambda)V_{-}(\omega,\lambda)-\mu_{-}(\lambda)V_{*}(\omega,\lambda)\right)}{\left(\mu_{*}(\lambda)V_{+}(\omega,\lambda)+\mu_{+}(\lambda)V_{*}(\omega,\lambda)\right)\left(\mu_{-}(\lambda)V_{*}(\omega,\lambda)+\mu_{*}(\lambda)V_{-}(\omega,\lambda)\right)}=:g(\lambda),$
(66)
$\begin{cases}&\mu_{*}(\lambda)V_{+}(\omega,\lambda)+\mu_{+}(\lambda)+V_{*}(\omega,\lambda)=0,\
\mu_{*}(\lambda)V_{-}(\omega,\lambda)-\mu_{-}(\lambda)V_{*}(\omega,\lambda)=0\\\
&\text{or}\\\
&\mu_{*}(\lambda)V_{+}(\omega,\lambda)-\mu_{+}(\lambda)V_{*}(\omega,\lambda)=0,\
\mu_{*}(\lambda)V_{-}(\omega,\lambda)+\mu_{-}(\lambda)V_{*}(\omega,\lambda)=0.\end{cases}$
(67)

Clearly, because $\omega\notin(M_{+}^{(k)}\cup M_{-}^{(k)}\cup
O^{(k)}\cup\Omega_{0})$, we get $\omega\notin(M_{+}^{(k,\lambda)}\cup
M_{-}^{(k,\lambda)}\cup O^{(k,\lambda)}\cup\Omega_{0})$ for all $\lambda$
small enough.

It remains to be shown that (66) cannot hold if $\lambda$ is nonzero and small
enough. Otherwise there would be a sequence $(\lambda_{j})\subset\mathbb{C}$
such that $\lambda_{j}\to 0$ and
$f(\lambda_{j}):=e^{2\mu_{*}(\lambda_{j})d}-g(\lambda_{j})=0\ \forall j.$ Due
to the $C^{1}$ nature of $f$, this would imply $f^{\prime}(0)=0.$ Because
$f^{\prime}(0)=2g(0)\mu_{*}^{\prime}(0)d-g^{\prime}(0)=-{\rm
i}g(0)d(\omega+V_{*}(\omega))\frac{1}{\mu_{*}(0)}-g^{\prime}(0)$, the equality
$f^{\prime}(0)=0$ is equivalent to

$d=\frac{{\rm i}\mu_{*}(0)g^{\prime}(0)}{g(0)(\omega+V_{*}(\omega))}.$

A direct calculation shows that $\frac{{\rm
i}\mu_{*}(0)g^{\prime}(0)}{g(0)(\omega+V_{*}(\omega))}=\beta(\omega)$. ∎

######  Example 1.

We choose the wave number $k=2$ and the $\mathcal{P}\mathcal{T}$-symmetric
setting (see Theorem 1.2 for the definition of
$\mathcal{P}\mathcal{T}$-symmetry) with three homogenous layers

$V_{*}(\omega)=-\frac{6}{5}\omega,\
V_{-}(\omega)=-\omega\left(1-\frac{2\pi}{\omega^{2}-{\rm
i}\frac{\omega}{2}}\right),\ V_{+}(\omega)=\overline{V_{-}(\omega)}.$

A width $d$ of the middle layer and corresponding eigenvalues can be found
using (40). In Figure 1 we plot $d_{j}(\omega)$ for $j=-1,0,1$ and
$\omega\in(0,5)$. Apparently $d_{-1}(\omega)\notin(0,\infty)$ for all
$\omega\in(0,5)$, $d_{0}(\omega)\in(0,\infty)$ for $\omega\in(0,a)\cup(b,5)$,
where $a\approx 2.5$, $b\approx 3.75$, and, finally,
$d_{1}(\omega)\in(0,\infty)$ for $\omega\in(c,5)$ with $c\approx 1.83$.

Figure 1: Functions $d_{j}(\omega)$ with $j=-1,0,1$ from (40).

As one can see, for example, for $d=0.7$ there are at least five eigenvalues
in $N^{(k)}$, see the intersections with the dotted line $d=0.7$ in Figure 1
(b) and (c). The approximate values are
$\omega_{0}\in\\{1.11,1.68,3.09,3.44,4.57\\}$, the first two of which are
visible in Figure 1 (b) and the last three in Figure 1 (c).

Next, we set the width of the middle layer equal to $d_{1}(\omega)$ and check
the algebraic simplicity and isolatedness of the corresponding eigenvalue
$\omega\in(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 $\alpha(\omega)$ and $|d_{1}(\omega)-\beta(\omega)|$ from
Propositions 2.7 and 2.8. Clearly, $\alpha$ remains positive on the whole
$(c,5)$ and $d_{1}\neq\beta$ everywhere except for at most two points, namely
$\omega\approx 2.15$ and $\omega\approx 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.

Figure 2: Quantities $\alpha(\omega)$ and $|d_{1}(\omega)-\beta(\omega)|$ with
$\alpha$ and $\beta$ from Propositions 2.7 and 2.8 respectively.

##  3 Bifurcation of nonlinear surface plasmons

The aim of this section is to prove Theorem 1.1. As stated in this theorem,
the solution of (5) is expanded in

$\omega=\omega_{0}+\varepsilon\nu+\varepsilon^{2}\sigma,\quad
u=\varepsilon^{1/2}\varphi_{0}+\varepsilon^{3/2}\phi+\varepsilon^{5/2}\psi,$
(68)

where $\omega_{0}\in\mathbb{C},\varphi_{0}\in D(A)$ are fixed by the linear
eigenvalue problem (and the normalization $\|\varphi_{0}\|=1$) and
$\nu\in\mathbb{C},\phi\in D(A)\cap\langle\varphi_{0}^{*}\rangle^{\perp}$ are
fixed by (13) and (14). Note that $\nu$ is well defined because of the
assumption (A-T) and that the normalization
$\langle\varphi_{0},\varphi_{0}^{*}\rangle=1$ is allowed because
$\langle\varphi_{0},\varphi_{0}^{*}\rangle=0$ would imply the solvability of
$L_{k}(\omega_{0})v=\varphi_{0}$ for $v\in D(A)$, which contradicts the
algebraic simplicity of $\omega_{0}$.

We decompose the nonlinear equation (5) using the projection

$P_{0}:D(A)\to\langle\varphi_{0}\rangle,u\mapsto\langle
u,\varphi_{0}^{*}\rangle\varphi_{0}$

and the complementary projection
$Q_{0}:=I-P_{0}:D(A)\to\langle\varphi_{0}^{*}\rangle^{\perp}$. This discussion
as well as the subsequent fixed point arguments are very similar to those in
[10]. We include them here (in a more compact form) for completeness.

Applying $P_{0}$ to equation (5), we get

$\begin{split}\langle L_{k}(\cdot,\omega)u,\varphi_{0}^{*}\rangle&=\langle
h(\cdot,\omega,u),\varphi_{0}^{*}\rangle\\\ &=\langle
h(\cdot,\omega,u)-\varepsilon^{\frac{3}{2}}h(\cdot,\omega_{0},\varphi_{0}),\varphi_{0}^{*}\rangle+\varepsilon^{\frac{3}{2}}\langle
h(\cdot,\omega_{0},\varphi_{0}),\varphi_{0}^{*}\rangle.\end{split}$ (69)

On the other hand, Taylor expanding $B(\omega)$ in $\omega_{0}$ up to order
two (using assumption (A-V)), we have

$\begin{split}\langle L_{k}(\cdot,\omega)u,\varphi_{0}^{*}\rangle&=\langle
L_{k}(\cdot,\omega_{0})u,\varphi_{0}^{*}\rangle+\langle\left(L_{k}(\cdot,\omega)-L_{k}(\cdot,\omega_{0})\right)u,\varphi_{0}^{*}\rangle\\\
&=\langle
u,L_{k}(\cdot,\omega_{0})^{*}\varphi_{0}^{*}\rangle-\langle\left(B(\cdot,\omega)-B(\cdot,\omega_{0})\right)u,\varphi_{0}^{*}\rangle\\\
&=-\left\langle\left(\partial_{\omega}B(\cdot,\omega_{0})(\omega-\omega_{0})+\frac{1}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})(\omega-\omega_{0})^{2}+I(\cdot,\omega)(\omega-\omega_{0})^{3}\right)u,\varphi_{0}^{*}\right\rangle,\end{split}$

where

$I(x,\omega):=\frac{1}{2\pi\mathrm{i}}\int_{\partial
B_{r}(\omega_{0})}\frac{B(x,z)}{(z-\omega_{0})^{3}(z-\omega)}\,{\rm d}z$

for any $r<\delta$ and $\omega\in B_{r}(\omega_{0})$. Inserting now the
expansions of $\omega$ and $u$ from (68) produces

$\begin{split}-\langle
L_{k}(\cdot,\omega)u,\varphi_{0}^{*}\rangle&=\varepsilon^{\frac{3}{2}}\nu\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle+\varepsilon^{\frac{5}{2}}\sigma\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0}),\varphi_{0}^{*}\rangle\\\
&\quad+\varepsilon^{\frac{5}{2}}\bigg{\langle}\bigg{(}\nu\partial_{\omega}B(\cdot,\omega_{0})\phi+\frac{\nu^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})\varphi_{0}\bigg{)},\varphi_{0}^{*}\bigg{\rangle}+\langle
v(\cdot,\nu,\sigma,\phi,\psi),\varphi_{0}^{*}\rangle,\end{split}$ (70)

where

$\begin{split}v(\cdot,\nu,\sigma,\phi,\psi):&=\varepsilon^{\frac{7}{2}}\left(\sigma\partial_{\omega}B(\cdot,\omega_{0})\phi+(\nu+\varepsilon\sigma)\partial_{\omega}B(\cdot,\omega_{0})\psi+\frac{\nu^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})(\phi+\varepsilon\psi)\right.\\\
&\left.\quad+\frac{1}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})(2\nu\sigma+\varepsilon\sigma^{2})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\right)+\varepsilon^{\frac{7}{2}}I(\cdot,\omega)(\nu+\varepsilon\sigma)^{3}(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\end{split}$
(71)

with $\omega=\omega_{0}+\varepsilon\nu+\varepsilon^{2}\sigma.$

The first inner product on the right hand side of (69) is
$o(\varepsilon^{\frac{3}{2}})$ due to the Lipschitz continuity of $h$, see
(77). Comparing now (69) and (70), the terms of order
$\varepsilon^{\frac{3}{2}}$ match if and only if we set

$\nu:=-\frac{\langle
h(\cdot,\omega_{0},\varphi_{0}),\varphi_{0}^{*}\rangle}{\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle},$
(72)

which is well-defined thanks to assumption (A-T). From the rest of (69)-(70)
we obtain

$\begin{split}\varepsilon^{\frac{5}{2}}\sigma\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle=&-\varepsilon^{\frac{5}{2}}\Big{\langle}\left(\nu\partial_{\omega}B(\cdot,\omega_{0})\phi+\frac{\nu^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})\varphi_{0}\right),\varphi_{0}^{*}\Big{\rangle}\\\
&-\langle v(\cdot,\nu,\sigma,\phi,\psi),\varphi_{0}^{*}\rangle-\langle
h(\cdot,\omega,u)-\varepsilon^{\frac{3}{2}}h\left(\cdot,\omega_{0},\varphi_{0}\right),\varphi_{0}^{*}\rangle.\end{split}$
(73)

Next, we apply $Q_{0}$ to (5). First, we write

$Q_{0}L_{k}(\cdot,\omega)u=Q_{0}L_{k}(\cdot,\omega_{0})Q_{0}(\varepsilon^{\frac{3}{2}}\phi+\varepsilon^{\frac{5}{2}}\psi)-Q_{0}\big{(}B(\cdot,\omega)-B(\cdot,\omega_{0})\big{)}u.$

The application of $Q_{0}$ to (5) produces

$\begin{split}&Q_{0}L_{k}(\cdot,\omega_{0})Q_{0}(\phi+\varepsilon\psi)=\varepsilon^{-\frac{3}{2}}\left(Q_{0}h(\cdot,\omega,u)+Q_{0}\big{(}B(\cdot,\omega)-B(\cdot,\omega_{0})\big{)}u\right)\\\
&=\varepsilon^{-\frac{3}{2}}Q_{0}h(\cdot,\omega,u)+Q_{0}\left((\nu+\varepsilon\sigma)\partial_{\omega}B(\cdot,\omega_{0})+\varepsilon\frac{(\nu+\varepsilon\sigma)^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})+\varepsilon^{2}I(\cdot,\omega)(\nu+\varepsilon\sigma)^{3}\right)(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi).\end{split}$
(74)

It remains to solve the system (73), (74) for $\sigma$ and $\psi$ with
$\phi\in D(A)\cap\langle\phi_{0}^{*}\rangle^{\perp}$ given by (14).

###  3.1 Proof of Theorem 1.1

We take advantage of the fact that the nonlinearity $h$ is well defined in
${\mathcal{H}}^{1}$ due to the algebra property of $H^{1}(a,b)$ for any
$(a,b)\subset\mathbb{R}$. In other words,
$h(\cdot,\omega,u)\in{\mathcal{H}}^{1}$ for any $u\in{\mathcal{H}}^{1}$, see
Lemma B.1.

First of all, we reduce equation (74) by eliminating the $O(1)$ part. Because
$h(x,\omega,\alpha u)=\alpha^{3}h(x,\omega,u)$ for all $\alpha\geq 0$ and
$(x,\omega,u)\in\mathbb{R}\times
B_{\delta}(\omega_{0})\times(D(A)\cap{\mathcal{H}}^{1})$, the term
$h(x,\omega,u)$ in (74) equals
$\varepsilon^{3/2}h(x,\omega,\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)$.
The $O(1)$ part of (74) thus holds if and only if

$Q_{0}L_{k}(\cdot,\omega_{0})Q_{0}\phi=Q_{0}(h(\cdot,\omega_{0},\varphi_{0})+\nu\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0}).$

Due to (72) we have
$P_{0}(h(\cdot,\omega_{0},\varphi_{0})+\nu\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0})=0$
and hence

$Q_{0}L_{k}(\cdot,\omega_{0})Q_{0}\phi=h(\cdot,\omega_{0},\varphi_{0})+\nu\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},$
(75)

which is equivalent to (14) in
$D(A)\cap{\mathcal{H}}^{1}\cap\langle\phi_{0}^{*}\rangle^{\perp}$ and uniquely
solvable by the following lemma. Note that the right hand side is in
${\mathcal{H}}^{1}$ because $\partial_{\omega}V(\cdot,\omega_{0})\in
W^{1,\infty}(I_{j})$ for all $j$ and $\varphi_{0}\in{\mathcal{H}}^{1}$ (and
hence by the algebra property, also
$h(\cdot,\omega,\varphi_{0})\in{\mathcal{H}}^{1}$). The regularity of
$\varphi_{0}$ follows from Lemma 2.2.

######  Lemma 3.1.

$(Q_{0}L_{k}Q_{0})^{-1}:Q_{0}{\mathcal{H}}^{1}\to Q_{0}({\mathcal{H}}^{1}\cap
D(A))$ is a bounded operator, i.e. $\|(Q_{0}L_{k}Q_{0})^{-1}\|\leq M_{L}$ with
$M_{L}>0$.

###### Proof.

The boundedness of the inverse $(Q_{0}L_{k}Q_{0})^{-1}:Q_{0}L^{2}\to
Q_{0}D(A)$ (with $D(A)$ equipped with the graph norm) follows from the closed
range theorem [13, Sec. VI.6]. Indeed, the operator $Q_{0}L_{k}Q_{0}$ is
Fredholm (see [19, Theorem IV.5.28], where the fact that $\lambda=0$ is a
simple isolated eigenvalue of the eigenvalue problem (15) is used).

The $H^{1}$ regularity on each layer holds for the second and the third
component of each element in $D(A)$ by definition of $D(A)$. The first
component of $u:=(Q_{0}L_{k}Q_{0})^{-1}r$ with $r\in Q_{0}{\mathcal{H}}^{1}$
is given by $u_{1}=\frac{{\rm i}}{V(\cdot,\omega)}(r_{1}-{\rm i}ku_{3})$,
which lies in ${\mathcal{H}}^{1}$ since $r_{1},u_{3}\in{\mathcal{H}}^{1}$ and
$1/V_{j}\in W^{1,\infty}(I_{j})$ for all $j$. ∎

######  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.

Having satisfied (75), the rest of (74) produces the following equation for
$(\sigma,\psi)\in\mathbb{C}\times(D(A)\cap{\mathcal{H}}^{1})$:

$\begin{split}Q_{0}L_{k}(\cdot,\omega_{0})Q_{0}\psi=&\varepsilon^{-1}Q_{0}\big{(}h(\cdot,\omega,\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)-h(x,\omega_{0},\varphi_{0})\big{)}+Q_{0}(\nu\partial_{\omega}B(\cdot,\omega_{0})(\phi+\varepsilon\psi))\\\
&+Q_{0}\left(\sigma\partial_{\omega}B(\cdot,\omega_{0})+\frac{(\nu+\varepsilon\sigma)^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})+\varepsilon
I(\cdot,\omega)(\nu+\varepsilon\sigma)^{3}\right)(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\\\
=:&R(\sigma,\psi),\end{split}$

where $\omega=\omega_{0}+\varepsilon\nu+\varepsilon^{2}\sigma$. We write this
as the fixed point problem

$\psi=G(\sigma,\psi):=(Q_{0}L_{k}(\omega_{0})Q_{0})^{-1}R(\sigma,\psi).$ (76)

The system (73), (76) is solved below for $\sigma$ and $\psi$ via a nested
Banach fixed point approach analogously to [10]. First we solve (76) for
$\psi\in B_{r_{2}}(0)\subset D(A)$ for each $\sigma\in
B_{r_{1}}(0)\subset\mathbb{C}$ fixed with $r_{1}>0$ arbitrary and a suitable
$r_{2}=r_{2}(r_{1})$. The radius $r_{2}$ is chosen so that
$G(\sigma,\cdot):B_{r_{2}}(0)\to B_{r_{2}}(0)$ is a contraction if
$\varepsilon>0$ is small enough.

Then, having obtained $\psi=\psi(\sigma)\in B_{r_{2}}(0)$, we shall solve
equation (73) for $\sigma$ in $B_{r_{1}}(0)\subset\mathbb{C}$ with a suitable
$r_{1}$. Note that the fixed point argument for equation (73) requires the
Lipschitz continuity of $\sigma\mapsto\psi(\sigma)$, which we verify below.

######  Remark 10.

The alternative approach of solving the system (73), (76) for
$(\sigma,\psi)\in B_{r_{1}}(0)\times
B_{r_{2}}(0)\subset\mathbb{C}\times(D(A)\cap{\mathcal{H}}^{1})$ via a fixed
point argument simultaneously is possible. However,
$R(\sigma_{1},\psi)-R(\sigma_{2},\psi)$ involves the term
$Q_{0}\partial_{\omega}B(\omega_{0})\varphi_{0}(\sigma_{1}-\sigma_{2})$, which
is $O(1)$ and prevents the contraction property. This can be avoided by
substituting for $\sigma_{1}-\sigma_{2}$ from (73) but this makes the approach
more technical.

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 $\varepsilon$, $r_{1}$, and $r_{2}$ for
$\varepsilon=\varepsilon(r_{1},r_{2})>0$ small enough.

Firstly, for all $|\sigma|\leq r_{1}$ and $\|\psi\|_{{\mathcal{H}}^{1}}\leq
r_{2}$, we have

$\displaystyle\|h(\cdot,\omega,\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)-h(\cdot,\omega_{0},\varphi_{0})\|_{{\mathcal{H}}^{1}}$
(77)
$\displaystyle\qquad\leq\|h(\cdot,\omega,\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)-h(\cdot,\omega_{,}\varphi_{0})\|_{{\mathcal{H}}^{1}}+\|h(\cdot,\omega,\varphi_{0})-h(\cdot,\omega_{0},\varphi_{0})\|_{{\mathcal{H}}^{1}}$
$\displaystyle\qquad\leq
c\varepsilon\left(\|\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi\|_{{\mathcal{H}}^{1}}^{2}\right)\|\phi+\varepsilon\psi\|_{{\mathcal{H}}^{1}}+L_{a}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}|\varepsilon\nu+\varepsilon^{2}\sigma|$
$\displaystyle\qquad\leq
c(\varepsilon+\varepsilon^{2}r_{2}+\varepsilon^{2}r_{1}).$

Here we have used Proposition B.1. Next, for each $j\in\\{1,\dots,m\\}$

$\displaystyle\|(\nu+\varepsilon\sigma)^{3}I(\cdot,\omega)\|_{W^{1,\infty}(I_{j})}$
$\displaystyle\leq|\nu+\varepsilon\sigma|^{3}\frac{\max\limits_{z\in\partial
B_{r}(\omega_{0})}\|B(\cdot,z)\|_{W^{1,\infty}(I_{j})}}{r^{2}(r-|\omega-\omega_{0}|)}\;\quad\text{for
all}\;\;\omega\in B_{r^{\prime}}(\omega_{0})\;\ \text{if}\;\
0<r^{\prime}<r<\delta$
$\displaystyle\leq|\nu+\varepsilon\sigma|^{3}\frac{2M+|\omega_{0}|+\frac{r}{2}}{r^{3}}\quad\text{for
all}\quad\omega\in B_{r/2}(\omega_{0}),$ $\displaystyle\leq C(1+\varepsilon
r_{1}),$ (78)

where

$M:=\max_{j\in\\{1,\dots,m\\}}\sup_{\omega\in
B_{r}(\omega_{0})}\|V(\cdot,\omega)\|_{W^{1,\infty}(I_{j})},$

$\|H\|_{W^{1,\infty}(I_{j})}$ for a matrix $H$ denotes the maximum of the
$W^{1,\infty}(I_{j})$-norm of all entries in $H$. We have assumed w.l.o.g.
that $M\geq|\omega_{0}|+|\delta|$ (in order to estimate the entry ${\rm
i}\omega$ in $B$ also by $M$).

Due to (A-V) and (78)

$\displaystyle\|\nu\partial_{\omega}B(\cdot,\omega_{0})(\phi+\varepsilon\psi)+\sigma\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\|_{{\mathcal{H}}^{1}}$
$\displaystyle\leq
c\left(|\sigma|\|\varphi_{0}\|_{{\mathcal{H}}^{1}}+\left(1+\varepsilon|\sigma|\right)\|\phi\|_{{\mathcal{H}}^{1}}+\varepsilon\left(1+|\sigma|\varepsilon\right)\|\psi\|_{{\mathcal{H}}^{1}}\right)$
$\displaystyle\leq
c(1+(1+\varepsilon)r_{1}+(\varepsilon+\varepsilon^{2}r_{1})r_{2})$
$\displaystyle\|(\nu+\varepsilon\sigma)^{2}\partial^{2}_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\|_{{\mathcal{H}}^{1}}$
$\displaystyle\leq c(1+\varepsilon r_{1}+\varepsilon^{2}r_{2})$
$\displaystyle\|\varepsilon
I(\cdot,\omega)(\nu+\varepsilon\sigma)^{3}(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi)\|_{{\mathcal{H}}^{1}}$
$\displaystyle\leq
c\varepsilon|\nu+\varepsilon\sigma|^{3}\left(\|\varphi_{0}\|_{{\mathcal{H}}^{1}}+\varepsilon\|\phi\|_{{\mathcal{H}}^{1}}+\varepsilon^{2}\|\psi\|_{{\mathcal{H}}^{1}}\right)$
$\displaystyle\leq c\varepsilon(1+\varepsilon r_{1}+\varepsilon^{2}r_{2}).$

As a result,

$\displaystyle\|R(\sigma,\psi)\|_{{\mathcal{H}}^{1}}\leq$ $\displaystyle
c_{0}\left(1+r_{1}+\varepsilon(r_{1}+r_{2})\right)\leq 2c_{0}(1+r_{1})$ (79)

for all $\sigma\in B_{r_{1}}(0)\subset\mathbb{C}$ and $\psi\in
B_{r_{2}}(0)\subset D(A)\cap{\mathcal{H}}^{1}$. The constant $c_{0}>0$ is
independent of $\varepsilon,r_{1},$ and $r_{2}$ if
$\varepsilon=\varepsilon(r_{1},r_{2})$ is small enough. The property
$G(\sigma,\cdot):B_{r_{2}}(0)\to B_{r_{2}}(0)$ is thus satisfied, e.g., with
$r_{2}=r_{2}(r_{1}):=2c_{0}M_{L}(1+r_{1})$ if $\varepsilon>0$ is small enough.

For the contraction property we similarly estimate

$\displaystyle\|R(\sigma,\psi_{1})-R(\sigma,\psi_{2})\|_{{\mathcal{H}}^{1}}\leq$
$\displaystyle
c\varepsilon^{-1}\left[\varepsilon^{2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}+\varepsilon^{4}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}^{2}+\varepsilon^{6}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}^{3}\right.$
$\displaystyle+\left.\varepsilon^{2}(1+\varepsilon
r_{1})\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}\right]$ $\displaystyle\leq$
$\displaystyle c\varepsilon\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}$

for all $\psi_{1},\psi_{2}\in B_{r_{2}}(0)$, $\sigma\in B_{r_{1}}(0)$, and for
$\varepsilon=\varepsilon(r_{1},r_{2})$ small enough.

Given $r_{1}>0$, we have thus found a unique $\psi=\psi(\sigma)\in
B_{r_{2}}(0)\subset D(A)\cap{\mathcal{H}}^{1}$ for any $\sigma\in
B_{r_{1}}(0),$ where $r_{2}=2c_{0}M_{L}(1+r_{1})$ and $\varepsilon>0$ is small
enough.

We proceed with the second step of the nested argument and solve (73) for
$\sigma$ after substituting $\psi=\psi(\sigma)$. For this we need the
Lipschitz continuity of $\sigma\mapsto\psi(\sigma)$ with respect to the
${\mathcal{H}}^{1}$ norm. From $\psi(\sigma)=G(\sigma,\psi(\sigma))$ and Lemma
3.1 we get for $\sigma_{1},\sigma_{2}\in B_{r_{1}}(0)$

$\|\psi(\sigma_{1})-\psi(\sigma_{2})\|_{{\mathcal{H}}^{1}}\leq
M_{L}\|R(\sigma_{1},\psi_{1})-R(\sigma_{2},\psi_{2})\|_{{\mathcal{H}}^{1}},$

where $\psi_{j}:=\psi(\sigma_{j}),j=1,2$. We denote also
$\omega_{j}:=\omega_{0}+\varepsilon\nu+\varepsilon^{2}\sigma_{j},$
$u_{j}:=\varepsilon^{1/2}\varphi_{0}+\varepsilon^{3/2}\phi+\varepsilon^{5/2}\psi_{j}$.

Analogously to [10] (see page 15), using (78), one estimates

$\|(\nu+\varepsilon\sigma_{1})^{3}I(\cdot,\omega_{1})-(\nu+\varepsilon\sigma_{2})^{3}I(\cdot,\omega_{2})\|_{W^{1,\infty}(I_{j})}\leq
c\varepsilon^{4}|\sigma_{1}-\sigma_{2}|\qquad\forall j\in\\{1,\dots,m\\},$
(80)

where $c$ has a cubic dependence on $r_{1}$.

Using again Proposition B.1 , we get

$\begin{split}\|h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{1})-h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{2})\|_{{\mathcal{H}}^{1}}&\leq
c\varepsilon^{2}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}},\\\
\|h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{2})-h(\cdot,\omega_{2},\varepsilon^{-1/2}u_{2})\|_{{\mathcal{H}}^{1}}&\leq
c\varepsilon^{2}L_{a}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}|\sigma_{1}-\sigma_{2}|\end{split}$
(81)

for all $\sigma_{1,2}\in B_{r_{1}}(0)$ and $\varepsilon>0$ small enough.
Hence, with the triangle inequalities

$\displaystyle\|h(\cdot,\omega_{1},\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{1})-h(\cdot,\omega_{2},\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{2})\|_{{\mathcal{H}}^{1}}$
$\displaystyle\leq\|h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{1})-h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{2})\|_{{\mathcal{H}}^{1}}+\|h(\cdot,\omega_{1},\varepsilon^{-1/2}u_{2})-h(\cdot,\omega_{2},\varepsilon^{-1/2}u_{2})\|_{{\mathcal{H}}^{1}}$

and

$\displaystyle\|\sigma_{1}\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{1})-\sigma_{2}\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{2})\|_{{\mathcal{H}}^{1}}$
$\displaystyle\qquad\leq\|\sigma_{1}\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{1})-\sigma_{1}\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{2})\|_{{\mathcal{H}}^{1}}$
$\displaystyle\qquad+\|(\sigma_{1}-\sigma_{2})\partial_{\omega}B(\cdot,\omega_{0})(\varphi_{0}+\varepsilon\phi+\varepsilon^{2}\psi_{2})\|_{{\mathcal{H}}^{1}},$

and similarly for the other terms in $R$, one obtains

$\displaystyle\|R(\sigma_{1},\psi_{1})-R(\sigma_{2},\psi_{2})\|_{{\mathcal{H}}^{1}}$
$\displaystyle\leq
c\left(\varepsilon\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}+\varepsilon
L_{a}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}|\sigma_{1}-\sigma_{2}|\right.$
$\displaystyle\quad\left.+\|\varphi_{0}\|_{{\mathcal{H}}^{1}}|\sigma_{1}-\sigma_{2}|+\varepsilon^{2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}+\varepsilon^{3}|\sigma_{1}-\sigma_{2}|\right).$

This clearly leads to

$\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}\leq
c\|\varphi_{0}\|_{{\mathcal{H}}^{1}}|\sigma_{1}-\sigma_{2}|$ (82)

for all $\sigma_{1,2}\in B_{r_{1}}(0)$ and $\varepsilon>0$ small enough with a
modified $c$. Equation 82 is the Lipschitz continuity of
$\sigma\mapsto\psi(\sigma)$ in ${\mathcal{H}}^{1}$.

Equation (73) is equivalent to

$\sigma=S(\sigma),$

where

$\displaystyle S(\sigma)$
$\displaystyle:=c_{B}\left[\Big{\langle}\left(\nu\partial_{\omega}B(\cdot,\omega_{0})\phi+\frac{\nu^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})\varphi_{0}\right),\varphi_{0}^{*}\Big{\rangle}+\varepsilon^{-\frac{5}{2}}\langle
v(\cdot,\nu,\sigma,\phi,\psi),\varphi_{0}^{*}\rangle\right.$
$\displaystyle\quad\left.+\varepsilon^{-\frac{5}{2}}\langle
h(\cdot,\omega,u)-\varepsilon^{\frac{3}{2}}h\left(\cdot,\omega_{0},\varphi_{0}\right),\varphi_{0}^{*}\rangle\right],$
$\displaystyle c_{B}$
$\displaystyle:=-\langle\partial_{\omega}B(\cdot,\omega_{0})\varphi_{0},\varphi_{0}^{*}\rangle^{-1}.$

Next, we find $r_{1}>0$ such that $S:B_{r_{1}}(0)\to B_{r_{1}}(0)$. Using
(A-V) and (78), we have for $\varepsilon$ small enough

$\displaystyle\|v(\cdot,\nu,\sigma,\phi,\psi)\|_{L^{2}(\mathbb{R})}$
$\displaystyle\leq c_{1}\varepsilon^{7/2},$
$\displaystyle\|h(\cdot,\omega,u)-\varepsilon^{\frac{3}{2}}h\left(\cdot,\omega_{0},\varphi_{0}\right)\|_{L^{2}(\mathbb{R})}$
$\displaystyle\leq
c_{2}\varepsilon^{5/2}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\phi+\varepsilon\psi\|_{{\mathcal{H}}^{1}}+c_{3}\varepsilon^{5/2}(|\nu|+\varepsilon|\sigma|)\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}+c_{4}\varepsilon^{7/2}$

with $c_{1},c_{2}$ dependent on $r_{1}$ but independent of $\varepsilon$, and
$c_{3}$ independent of $r_{1}$ and $\varepsilon$ if
$\varepsilon=\varepsilon(r_{1})>0$ is small enough. Because
$\|\psi\|_{{\mathcal{H}}^{1}}\leq r_{2}=2c_{0}M_{L}(1+r_{1})$, this leads to

$|S(\sigma)|\leq
2c_{B}\left[c_{5}+c_{2}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\phi\|_{{\mathcal{H}}^{1}}+c_{3}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}(|\nu|+\varepsilon
r_{1})\right]$

for all $\varepsilon$ small enough, where
$c_{5}:=2\left|\Big{\langle}\left(\nu\partial_{\omega}B(\cdot,\omega_{0})\phi+\frac{\nu^{2}}{2}\partial_{\omega}^{2}B(\cdot,\omega_{0})\varphi_{0}\right),\varphi_{0}^{*}\Big{\rangle}\right|$.
Hence, we can choose, e.g.,
$r_{1}:=2c_{B}(c_{5}+c_{2}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\phi\|_{{\mathcal{H}}^{1}}+c_{3}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}|\nu|+1)$.

Let us now check the contraction property of $S$ on
$B_{r_{1}}(0)\subset\mathbb{C}$. Using (80) and (A-V), one easily estimates

$\displaystyle|\langle
v(\cdot,\nu,\sigma_{1},\phi,\psi_{1})-v(\cdot,\nu,\sigma_{2},\phi,\psi_{2}),\varphi_{0}^{*}\rangle|$
$\displaystyle\qquad\leq
c\left(\varepsilon^{7/2}|\sigma_{1}-\sigma_{2}|+\varepsilon^{9/2}|\sigma_{1}^{2}-\sigma_{2}^{2}|+\varepsilon^{13/2}|\sigma_{1}^{3}-\sigma_{2}^{3}|+\varepsilon^{7/2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}\right).$

From (81) we obtain

$|\langle
h(\cdot,\omega_{1},u_{1})-h(\cdot,\omega_{2},u_{2}),\varphi_{0}^{*}\rangle|\leq
c\varepsilon^{7/2}\left(\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{2}\|\psi_{1}-\psi_{2}\|_{{\mathcal{H}}^{1}}+L_{a}\|\varphi_{0}\|_{{\mathcal{H}}^{1}}^{3}|\sigma_{1}-\sigma_{2}|\right)$

with the above definitions of $\omega_{1},\omega_{2}$. Together with (82) the
last two inequalities lead to

$|S(\sigma_{1})-S(\sigma_{2})|\leq c\varepsilon|\sigma_{1}-\sigma_{2}|$

for all $\sigma_{1},\sigma_{2}\in B_{r_{1}}(0)$ and $\varepsilon$ small
enough. The contraction property follows for $\varepsilon>0$ small enough and
the proof is finished.

##  4 Proof of Theorem 1.2 (bifurcation under the
$\mathcal{P}\mathcal{T}$-symmetry)

We prove next Theorem 1.2, i.e. we show that the under the assumption of
$\mathcal{P}\mathcal{T}$-symmetry of the coefficients and the realness of
$\omega_{0}$, the above bifurcation argument can be carried out in the
$\mathcal{P}\mathcal{T}$-symmetric subspace resulting in real $\omega$ and
$\mathcal{P}\mathcal{T}$-symmetric $u$.

As we show now, after multiplication with $-{\rm i}$ the operator
$L_{k}(\cdot,\omega)$ as well as the nonlinearity $h$ commute with the
$\mathcal{P}\mathcal{T}$ symmetry provided $\hat{\chi}^{(1,3)}$ are
$\mathcal{P}\mathcal{T}$-symmetric. We have

$\tilde{L}_{k}(\cdot,\omega)u=\tilde{h}(\cdot,\omega,u(\cdot))$ (83)

for $u=\left(u_{1},u_{2},u_{3}\right)^{\mathsf{T}}$, where
$\tilde{L}_{k}(x,\omega):=-{\rm i}L_{k}(x,\omega)$ and
$\tilde{h}(x,\omega,u):=-{\rm i}h(x,\omega,u)$. For reference, we recall

$\begin{split}\tilde{L}_{k}(x,\omega)&=\begin{pmatrix}0\ &0\ &k\\\ 0\ &0\
&{\rm i}\partial_{x}\\\ -k\ &-\mathrm{i}\partial_{x}\
&0\end{pmatrix}-\begin{pmatrix}V(x,\omega)\ &0\ &0\\\ 0\ &V(x,\omega)\ &0\\\
0\ &0\ &\omega\end{pmatrix},\quad
V(x,\omega)=-\omega\varepsilon_{0}\mu_{0}\left(1+\hat{\chi}^{(1)}(x,\omega)\right),\\\
\tilde{h}(x,\omega,u)&=-\omega\epsilon_{0}\mu_{0}^{3}\hat{\chi}^{(3)}(x,\omega)\left(2|u_{E}|^{2}u_{E}+(u_{E}\cdot
u_{E})\overline{u_{E}}\right)\quad\text{with}\quad
u_{E}=\left(u_{1},u_{2},0\right)^{\mathsf{T}}.\end{split}$ (84)

######  Lemma 4.1.

Let $\omega\in\mathbb{R}$. If the functions $\hat{\chi}^{(1)}(x,\omega)$ and
$\hat{\chi}^{(3)}(x,\omega)$ are $\mathcal{P}\mathcal{T}$-symmetric, i.e.

$\hat{\chi}^{(1)}(x,\omega)=\overline{\hat{\chi}^{(1)}(-x,\omega)},\
\hat{\chi}^{(3)}(x,\omega)=\overline{\hat{\chi}^{(3)}(-x,\omega)}\ \forall
x\in\mathbb{R},$

then the operator $\tilde{L}_{k}(\omega)$ and the nonlinearity
$\tilde{h}(\cdot,\omega,u)$ commute with the
$\mathcal{P}\mathcal{T}$-symmetry, i.e.

$\mathcal{B}\tilde{L}_{k}(x,\omega)=\tilde{L}_{k}(x,\omega)\mathcal{B},\quad\mathcal{B}\tilde{h}(x,\omega,u(x))=\tilde{h}(x,\omega,\mathcal{B}u(x))\quad\text{with}\quad\mathcal{B}u(x):=\overline{u(-x)}.$

###### Proof.

Let $u=(u_{1},u_{2},u_{3})^{T}\in D(A)$. Due to the symmetry of $V(x,\omega)$
we get

$\begin{split}\mathcal{B}\tilde{L}_{k}(x,\omega)u&=\begin{pmatrix}k\overline{u_{3}(-x)}-\overline{V(-x,\omega)}\
\overline{u_{1}(-x)}\\\
-\mathrm{i}(\partial_{x}\overline{u_{3}})(-x)-\overline{V(-x,\omega)}\
\overline{u_{2}(-x)}\\\ -k\overline{u_{1}(-x)}+{\rm
i}(\partial_{x}\overline{u_{2}})(-x)-\omega\overline{u_{3}(-x)}\end{pmatrix}=\begin{pmatrix}k\overline{u_{3}(-x)}-V(x,\omega)\
\overline{u_{1}(-x)}\\\
\mathrm{i}\partial_{x}(\overline{u_{3}(-x)})-V(x,\omega)\
\overline{u_{2}(-x)}\\\ -k\overline{u_{1}(-x)}-{\rm
i}\partial_{x}(\overline{u_{2}(-x)})-\omega\overline{u_{3}(-x)}\end{pmatrix}\\\
&=\tilde{L}_{k}(x,\omega)\mathcal{B}u.\end{split}$

Hence $\mathcal{B}\tilde{L}_{k}(\omega)=\tilde{L}_{k}(\omega)\mathcal{B}$.

Next, we show the $\mathcal{P}\mathcal{T}$-symmetry of the nonlinearity
$\tilde{h}$. As $\hat{\chi}^{(3)}(x,\omega)$ is
$\mathcal{P}\mathcal{T}$-symmetric, one, indeed, obtains

$\begin{split}\mathcal{B}(\tilde{h}(\cdot,\omega,u))(x)&=-\omega\varepsilon_{0}\mu_{0}^{3}\overline{\hat{\chi}^{(3)}(-x,\omega)}\
\left(\overline{2|u_{E}|^{2}u_{E}+(u_{E}\cdot
u_{E})\overline{u_{E}}}\right)(-x)\\\
&=-\omega\varepsilon_{0}\mu_{0}^{3}\hat{\chi}^{(3)}(x,\omega)\
\left(2|u_{E}(-x)|^{2}\overline{u_{E}(-x)}+(\overline{u_{E}(-x)}\cdot\overline{u_{E}(-x)})u_{E}(-x)\right)\\\
&=\tilde{h}(\cdot,\omega,\mathcal{B}u)(x).\end{split}$

∎

Lemma 4.1 allows us to restrict the bifurcation problem to the
$\mathcal{P}\mathcal{T}$-symmetric subspace of $D(A)\cap{\mathcal{H}}_{1}$,
i.e. to $\\{u\in D(A)\cap{\mathcal{H}}_{1}:\mathcal{B}u=u\\}$ provided the
$\mathcal{P}\mathcal{T}$-symmetry of $\hat{\chi}^{(1,3)}(\cdot,\omega)$ holds
for all $\omega$ in a real neighbourhood of $\omega_{0}$. The bifurcating
eigenvalue is then real and the solution $u$ is
$\mathcal{P}\mathcal{T}$-symmetric, i.e. Theorem 1.2 holds. This is proved
completely analogously to Proposition 3.1 in [10].

Note that the proof uses the $\mathcal{P}\mathcal{T}$-symmetry of the
eigenfunction $\varphi_{0}$. Under the assumption of algebraic simplicity and
realness of $\omega_{0}$ the eigenfunction can always be selected
$\mathcal{P}\mathcal{T}$-symmetric. Indeed, applying $\mathcal{B}$ to the
eigenvalue equation yields $L_{k}(\omega_{0})(\mathcal{B}\varphi_{0})=0$ due
to the $\mathcal{P}\mathcal{T}$-symmetry of $L_{k}$. As $\omega_{0}$ is
simple, we get $\varphi_{0}=\mathcal{B}\varphi_{0}$ (up to a multiplicative
constant).

##  5 Numerical Results

We restrict the numerical computations to two layers ($m=2$) with the
interface at $x=0$. We first reduce the nonlinear system (5) with $\omega\neq
0$ to a system of two equations for $\tilde{u}:=(u_{1},u_{2})$. For
$\omega\neq 0$ the third equation in (5) namely produces

$u_{3}=-\frac{{\rm i}}{\omega}\left(u_{2}^{\prime}-{\rm i}ku_{1}\right)$ (85)

and the remaining equations become

$\displaystyle\tilde{u}_{2}^{\prime}-{\rm i}\left(k+\frac{\omega
V(x)}{k}\right)\tilde{u}_{1}-\frac{\omega}{k}h_{1}(x,\omega,\tilde{u})$
$\displaystyle=0,\quad x\in\mathbb{R}\setminus\\{0\\},$ (86)
$\displaystyle\tilde{u}_{2}^{\prime\prime}-{\rm i}k\tilde{u}_{1}-\omega
V(x)\tilde{u}_{2}+{\rm i}\omega h_{2}(x,\omega,\tilde{u})$
$\displaystyle=0,\quad x\in\mathbb{R}\setminus\\{0\\}.$ (87)

Due to equation (85) the interface conditions (11) become

$\llbracket\tilde{u}_{2}\rrbracket=\llbracket{\rm
i}k\tilde{u}_{1}-\tilde{u}_{2}^{\prime}\rrbracket=0\ \text{ at }x=0.$ (88)

As (86)-(88) is to be solved for $(\tilde{u}_{1},\tilde{u}_{2})$ and $\omega$,
an additional constraint is needed. We choose the condition

$\langle u,\varphi_{0}^{*}\rangle=\sqrt{\varepsilon}$ (89)

with a given $\varepsilon>0$. By varying $\varepsilon$ in a right
neighbourhood of zero we will thus generate a bifurcation curve in accordance
with Theorem 1.1. System (86)-(89) is solved via the Newton iteration in a
finite difference discretization. Because the function
$|\tilde{u}|^{2}\tilde{u}$ is not complex differentiable, we work in the real
variables

$\tilde{u}_{1,R},\tilde{u}_{2,R},\tilde{u}_{1,I},\tilde{u}_{2,I},\omega_{R},\
\text{and}\ \omega_{I},$

where
$\tilde{u}_{j,R}=\operatorname{Re}(\tilde{u}_{j}),\tilde{u}_{j,I}=\operatorname{Im}(\tilde{u}_{j}),j=1,2$,
$\omega_{R}=\operatorname{Re}(\omega)$, and
$\omega_{I}=\operatorname{Im}(\omega)$. We rewrite (86)-(89) as a system of
six real equations with four real interface conditions, namely as the real and
imaginary parts of equations (86)-(89).

###  5.1 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=\pm L$, which are suitable as we search for localized solutions.

Next, we explain that while equation (87) is to be solved at all grid points
excluding $x=0$, equation (86) has to be solved also in the limit $x\to 0+$
and $x\to 0-$. The reason is the condition $\nabla\cdot D=0$. The divergence
condition

$0=\nabla\cdot D=-\frac{1}{\omega}\left[\partial_{x}\left(V\tilde{u}_{1}-{\rm
i}h_{1}(x,\omega,\tilde{u})\right)+{\rm i}k\left(V\tilde{u}_{1}-{\rm
i}h_{2}(x,\omega,\tilde{u})\right)\right]$

is satisfied on $\mathbb{R}\setminus\\{0\\}$ by solutions of (86), (87)
because these equations imply

$V\tilde{u}_{1}-{\rm i}h_{1}=\frac{{\rm i}k}{\omega}\left({\rm
i}k\tilde{u}_{1}-\tilde{u}_{2}^{\prime}\right),\quad V\tilde{u}_{2}-{\rm
i}h_{2}=\frac{1}{\omega}\left(\tilde{u}_{2}^{\prime\prime}-{\rm
i}k\tilde{u}_{1}\right).$

To get $\nabla\cdot D=0$ distributionally on $\mathbb{R}$ (which is the
correct formulation of the divergence condition for $u\in D(A)$), we need to
satisfy also $\llbracket D_{1}\rrbracket=0,$ i.e. $\llbracket
V\tilde{u}_{1}-{\rm i}h_{1}\rrbracket=0$. This follows automatically from the
limits $x\to 0\pm$ of (86) and the second interface condition in (88). After
discretization it means that we need to solve (86) also in the limits $x\to
0-$ and $x\to 0+$.

Choosing $\Delta x:=\frac{2L}{N+1}$ with $N\in 2\mathbb{N}+1$, we have the $N$
grid points $-\frac{N-1}{2}\Delta x,-\frac{N-2}{2}\Delta x,\dots-\Delta
x,0,\Delta x,\dots,\frac{N-1}{2}\Delta x$. We denote
$x_{j}:=-\frac{N-1}{2}\Delta x+(j-1)\Delta x$, $j=1,\dots,N$. Let
$j_{*}:=\frac{N+1}{2}$ be the index of the interface grid point, i.e.
$x_{j_{*}}=0$.

Because $\tilde{u}_{1}$ is discontinuous at $x=0$, the degrees of freedom of
the discretized problem must include approximations of $\tilde{u}_{1}(0-)$ and
$\tilde{u}_{1}(0+)$. Hence, we have the (real) degrees of freedom

$U^{R/I}_{1,j},U^{R/I}_{2,j},\quad x\in\\{1,\dots,j_{*}-1,j_{*}+1,\dots,N\\}$

approximating $\tilde{u}_{1,R/I}$ at $x_{j},j\neq j_{*}$ and

$U^{R/I}_{1,j_{*},-},U^{R/I}_{1,j_{*},+},U^{R/I}_{2,j_{*}},$

where $U^{R/I}_{1,j_{*},\pm}$ approximate $\tilde{u}_{1,R/I}(0\pm)$ and
$U^{R/I}_{2,j_{*}}$ approximate $\tilde{u}_{2,R/I}(0)$. Finally, there are the
two degrees of freedom $\omega_{R}$ and $\omega_{I}$ approximating
$\operatorname{Re}(\omega)$ and $\operatorname{Im}(\omega)$.

These are altogether $4N+4$ real degrees of freedom. Solving (86) in the
finite difference discretization at $x_{j},j\neq j_{*}$ as well as at
$x_{1}\to 0\pm$ and (87) at $x_{j},j\neq j_{*}$, we get $2N+2$ real equations
from (86) and $2N-2$ equations from (87). Condition (89) produces two real
equations. After using the interface condition $\llbracket{\rm
i}k\tilde{u}_{1}-\tilde{u}_{2}^{\prime}\rrbracket=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 $x_{j},j\neq j_{*}$, we use the centered finite difference stencil of
second order for $\partial_{x}$ as well as for $\partial_{x}^{2}$. At $x=0-$
and $x=0+$ we use the one-sided second order stencils. The integral in the
normalization condition (89) is approximated using the trapezoidal rule.

###  5.2 Bifurcation Examples

We present here two numerical examples of bifurcation in the case of 2 layers.
One of the examples is $\mathcal{P}\mathcal{T}$-symmetric and the linear
frequency $\omega_{0}$ (hence also the bifurcating frequency $\omega$) 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
$\omega=\omega_{0}+\varepsilon\nu$ with $0<\varepsilon\ll 1$ is provided by
the asymptotic approximation $\varepsilon^{1/2}\varphi_{0}$, see Theorem 1.1.
The bifurcation curve is obtained by a simple parameter continuation in
$\varepsilon$.

####  5.2.1 $\mathcal{P}\mathcal{T}$-symmetric example

We choose $k=1$ and a case of a $\mathcal{P}\mathcal{T}$-symmetric Drude
material that is homogenous on each layer $x<0$ and $x>0$. In detail,

$V_{-}(\omega)=-\omega\left(1-\frac{2\pi\omega_{p}^{2}}{\omega^{2}+{\rm
i}\gamma\omega}\right),\quad V_{+}=\overline{V_{-}}$ (90)

with $\omega_{p}=0.5$ and $\gamma=0.7$. The resulting function $V$ is plotted
in Figure 3 (a).

We find a real eigenvalue $\omega_{0}$ by determining real elements of
$N^{(k)}$. In the case of (90) condition (21) can, in fact, be solved
explicitly. First, we note that (21) is equivalent to
$k^{2}(W_{+}(\omega)+W_{-}(\omega))=W_{+}(\omega)W_{-}(\omega)$. This
equation, with $W\pm=-\omega V_{\pm}(\omega)$ and $V_{\pm}$ given in (90),
reduces to a quadratic equation in $\mu:=\omega^{2}$. The two roots are

$\mu_{1,2}:=k^{2}+2\pi\omega_{p}^{2}-\frac{\gamma^{2}}{2}\pm\frac{1}{2}\left(\gamma^{2}(\gamma^{2}-4(2\pi\omega_{p}^{2}-k^{2}))+4k^{4}\right)^{1/2}$

and we get $\mu_{1}=\omega_{0}^{2}$ with $\omega_{0}\approx 1.791$. The
corresponding eigenfunction $\varphi_{0}$ is plotted in Figure 4.

In the nonlinearity $h$ we choose $\epsilon_{0}\mu_{0}^{3}\chi^{(3)}\equiv 1,$
i.e. $\chi^{(3)}$ is real and $x-$ as well as $\omega-$independent. In
particular, it is also $\mathcal{P}\mathcal{T}$ symmetric.

The discretization parameters are $L=120$ and $N=17999$ (resulting in $\Delta
x\approx 0.013$). The numerically approximated value of $\nu$ is
$\approx-0.2572$.

In Figure 3 (b) we plot the bifurcation diagram showing that the bifurcation
parameter (i.e. frequency) $\omega$ indeed stays real. As expected, the
bifurcation curve is tangent to the line given by the asymptotic approximation
$\omega=\omega_{0}+\nu\varepsilon$ with $\nu$ from (13). In Figure 3 (c) we
show that the convergence of the asymptotic approximation error
$|\omega-\omega_{0}-\varepsilon\nu|$ is indeed approximately quadratic. The
solution $u$ (chosen at $\omega\approx 1.7167$) plotted in 4 satisfies the
$\mathcal{P}\mathcal{T}$-symmetry. It is clearly close to the linear solution
$\varphi_{0}$ but not identical.

Figure 3: (a) The graph of the $\mathcal{P}\mathcal{T}$-symmetric
$V(\cdot,\omega_{0})$ given by (90) with $\omega=\omega_{0}\approx 1.7914$.
(b) The corresponding bifurcation diagram for the bifurcation from the
eigenvalue $\omega_{0}$. (c) Convergence of the approximation error
$|\omega-\omega_{0}-\varepsilon\nu|$. For comparison a curve with a quadratic
convergence is plotted. Figure 4: The nonlinear solution $u$ at $\omega\approx
1.7167$. Recall that $E_{1}=u_{1},E_{2}=u_{2},H_{3}=u_{3}$. The linear
eigenfunction $\varphi_{0}$ (normalized ot have a similar amplitude to that of
$u$) is plotted for comparison.

####  5.2.2 Non-$\mathcal{P}\mathcal{T}$-symmetric example

We present also one example which is not $\mathcal{P}\mathcal{T}$-symmetric,
namely

$V_{-}(\omega)=-\omega\left(1-\frac{2\pi\omega_{p}^{2}}{\omega^{2}+{\rm
i}\gamma\omega}\right),\quad V_{+}=-\omega(1+\eta)$ (91)

with $\omega_{p}=0.8,\gamma=1$, and $\eta=1$. Because of the lack of
$\mathcal{P}\mathcal{T}$-symmetry, the linear eigenvalues are not real.
Equation (5) makes sense also for non-real omega. 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
$\omega_{0}\in\mathbb{C}\setminus\mathbb{R}$.

In the nonlinearity $h$ we choose again
$\epsilon_{0}\mu_{0}^{3}\chi^{(3)}\equiv 1.$ The discretization parameters are
$L=100$ and $N=11999$. The value of $\nu$ is approximated as
$\nu\approx-0.0336-{\rm i}0.0054.$

Solving (21) numerically, we obtain an eigenvalue $\omega_{0}\approx
0.4679-{\rm i}\leavevmode\nobreak\ 0.061.$ In Fig. 5 (a) we plot $V$ given by
(91) with $\omega=\omega_{0}$. In (b) we show the bifurcation diagram and the
first order asymptotic approximation $\omega=\omega_{0}+\nu\varepsilon$. The
corresponding asymptotic error is plotted in (c) with an observed
approximately quadratic convergence, as predicted.

Figure 5: (a) The graph of the $V(\cdot,\omega_{0})$ given by (91) with
$\omega=\omega_{0}\approx 0.4679-{\rm i}\leavevmode\nobreak\ 0.061$. (b) The
corresponding bifurcation diagram for the bifurcation from the eigenvalue
$\omega_{0}$. (c) Convergence of the approximation error
$|\omega-\omega_{0}-\varepsilon\nu|$. For comparison a curve with a quadratic
convergence is plotted.

The nonlinear solution $u$ at $\omega\approx 0.427-{\rm i}\leavevmode\nobreak\
0.066$ is plotted in (b) together with the linear eigenfunction $\varphi_{0}$
normalized to have a similar amplitude to that of $u$.

Figure 6: The nonlinear solution $u$ at $\omega\approx 0.4276-{\rm i}0.066$.
Recall that $E_{1}=u_{1},E_{2}=u_{2},H_{3}=u_{3}$. The linear eigenfunction
$\varphi_{0}$ (normalized ot have a similar amplitude to that of $u$) is
plotted for comparison.

##  Appendix A Solvability of $L_{k}(\omega_{0})u=r$ with
$\omega_{0}\in\sigma_{p}(L_{k})$ and $r\in Q_{0}L^{2}(\mathbb{R})$ for the
case of two homogenous layers using the variation of parameters

Here we consider the linear case of two homogenous layers of Section 2.1.1,
i.e. $m=2$, $V_{1}(x,\omega)=V_{-}(\omega),V_{2}(x,\omega)=V_{+}(\omega)$ with
the interface at $x=0$. Let $\omega_{0}\in\sigma_{p}(L_{k})$ be simple and
$\varphi_{0}\in D(A)$, $\varphi_{0}^{*}\in D(A^{*})=D(A)$ be corresponding
eigenfunctions, i.e.
$L_{k}(\cdot,\omega_{0})\varphi_{0}=0,L_{k}^{*}(\cdot,\omega_{0})\varphi_{0}^{*}=0$.
As we know, using the closed range theorem, equation

$L_{k}(x,\omega)u=r\quad\text{ with }\quad\langle r,\varphi_{0}^{*}\rangle=0$

has a unique solution in $D(A)\cap\langle\varphi_{0}^{*}\rangle^{\perp}$. Here
we want to demonstrate this explicitly using the variation of parameters.

The solution $u\in D(A)$ is given by (27), (29) and

$\left(\mu_{-}V_{+}(\omega_{0})+\mu_{+}V_{-}(\omega_{0})\right)C_{-}=\frac{1}{2}\left(\mu_{+}V_{-}(\omega_{0})-\mu_{-}V_{+}(\omega_{0})\right)\int_{-\infty}^{0}\rho_{-}^{(2)}(s)e^{\mu_{-}s}\,{\rm
d}s+\mu_{+}V_{+}(\omega_{0})\int_{0}^{\infty}\rho_{+}^{(1)}(s)e^{-\mu_{+}s}\,{\rm
d}s$ (92)

with $\rho_{+}^{(1)}$ and $\rho_{-}^{(2)}$ in (26).

Because $\omega_{0}\in N^{(k)}$, the left hand side of (92) vanishes, see
(21). The existence of $u\in D(A)$ follows if the right hand side vanishes
too. This is shown below to hold if $r\in Q_{0}L^{2}(\mathbb{R})$, i.e. if
$\langle r,\varphi_{0}^{*}\rangle=0.$

From

$L_{k}^{*}(\cdot,\omega_{0})=A-\overline{B(\cdot,\omega_{0})}=\overline{L_{-k}(\cdot,\omega_{0})},$

we conclude

$\varphi_{0}^{*}(x;k)=\overline{\varphi_{0}(x;-k)}.$

The eigenfunction $\varphi_{0}(x;k)$ is given by (35) with
$c_{+}=-\frac{V_{-}(\omega_{0})}{V_{+}(\omega_{0})}c_{-},$
$\mu_{\pm}=\sqrt{k^{2}-W_{\pm}(\omega_{0})}$, where
$V_{+}(\omega_{0})\mu_{-}+V_{-}(\omega_{0})\mu_{+}=0.$ Hence,

$\varphi_{0}^{*}(x;k)=\begin{cases}\overline{c_{-}}e^{\overline{\mu_{-}}x}\left(\begin{smallmatrix}-{\rm
i}k\\\ \overline{\mu_{-}}\\\ {\rm
i}\overline{V_{-}(\omega)}\end{smallmatrix}\right)&\quad\text{for}\ x<0,\\\
\overline{c_{+}}e^{-\overline{\mu_{+}}x}\left(\begin{smallmatrix}{\rm i}k\\\
\overline{\mu_{+}}\\\ -{\rm
i}\overline{V_{+}(\omega)}\end{smallmatrix}\right)&\quad\text{for}\
x>0.\end{cases}$

The assumption $\langle r,\varphi_{0}^{*}\rangle=0$ becomes

$\displaystyle\int_{-\infty}^{0}e^{\mu_{-}x}\left({\rm
i}kr_{1}(x)+\mu_{-}r_{2}(x)-{\rm i}V_{-}(\omega_{0})r_{3}(x)\right)\,{\rm d}x$
(93)
$\displaystyle-\frac{V_{-}(\omega_{0})}{V_{+}(\omega_{0})}\int_{0}^{\infty}e^{-\mu_{+}x}\left(-{\rm
i}kr_{1}(x)+\mu_{+}r_{2}(x)+{\rm i}V_{+}(\omega_{0})r_{3}(x)\right)\,{\rm
d}x=0.$

The right hand side of (92) equals

$\displaystyle\mu_{+}V_{-}(\omega_{0})$
$\displaystyle\int_{-\infty}^{0}\left(\frac{-kr_{1}(x)}{V_{-}(\omega_{0})\mu_{-}}+{\rm
i}\frac{r_{2}(x)}{V_{-}(\omega_{0})}+\frac{r_{3}(x)}{\mu_{-}}\right)e^{\mu_{-}x}\,{\rm
d}x$
$\displaystyle+\mu_{+}V_{+}(\omega_{0})\int_{0}^{\infty}\left(\frac{kr_{1}(x)}{V_{+}(\omega_{0})\mu_{+}}+{\rm
i}\frac{r_{2}(x)}{V_{+}(\omega_{0})}-\frac{r_{3}(x)}{\mu_{+}}\right)e^{-\mu_{+}x}\,{\rm
d}x$

and using (93), it simplifies to

$\left(\frac{\mu_{+}V_{-}(\omega_{0})}{\mu_{-}V_{+}(\omega_{0})}+1\right)\int_{0}^{\infty}\left(\frac{kr_{1}(x)}{V_{+}(\omega_{0})\mu_{+}}+{\rm
i}\frac{r_{2}(x)}{V_{+}(\omega_{0})}-\frac{r_{3}(x)}{\mu_{+}}\right)e^{-\mu_{+}x}\,{\rm
d}x,$

which is indeed zero since
$V_{+}(\omega_{0})\mu_{-}+V_{-}(\omega_{0})\mu_{+}=0.$

##  Appendix B Lipschitz continuity of the nonlinearity $h$

In the following we use the notation $\lesssim$ to denote the inequality up to
a multiplicative constant, which is independent of the variables and functions
appearing on the right hand side of the estimate.

######  Proposition B.1.

Let $\delta>0$ be such that $B_{\delta}(\omega_{0})\subset\Omega$ and such
that (A-Na) holds. Then the nonlinearity $h$ given by

$h(x,\omega,u)=-{\rm
i}\varepsilon_{0}\mu_{0}^{3}\hat{\chi}^{(3)}(x,\omega)\left(2|\tilde{u}|^{2}(\tilde{u},0)^{\mathsf{T}}+(\tilde{u}\cdot\tilde{u})\left(\overline{\tilde{u}},0\right)^{\mathsf{T}}\right),\quad\text{where}\quad
u=(\tilde{u},u_{3})^{\mathsf{T}},\ \tilde{u}=(u_{1},u_{2}),$

satisfies $h(\cdot,\omega,u)\in{\mathcal{H}}^{1}$ for all $u\in
D(A)\cap{\mathcal{H}}^{1}$ and $\omega\in B_{\delta}(\omega_{0})$ as well as
the Lipschitz properties

$\left\|h(\cdot,\omega,\psi)-h(\cdot,\omega,\varphi)\right\|_{{\mathcal{H}}^{1}}\lesssim\left(\left\|\psi\right\|^{2}_{{\mathcal{H}}^{1}}+\left\|\varphi\right\|^{2}_{{\mathcal{H}}^{1}}\right)\left\|\psi-\varphi\right\|_{{\mathcal{H}}^{1}},$
(94)

for all $\omega\in B_{\delta}(\omega_{0})$ and $\varphi$, $\psi\in
D(A)\cap{\mathcal{H}}^{1}$;

$\left\|h(\cdot,\omega_{1},u)-h(\cdot,\omega_{2},u)\right\|_{{\mathcal{H}}^{1}}\lesssim\left\|u\right\|^{3}_{{\mathcal{H}}^{1}}\left|\omega_{1}-\omega_{2}\right|$
(95)

for all $\omega_{1},\omega_{2}\in B_{\delta}(\omega_{0})$ and $u\in
D(A)\cap{\mathcal{H}}^{1}$.

###### Proof.

Due to the algebra property of ${\mathcal{H}}^{1}$ it follows that

$\begin{split}\left\|h(\cdot,\omega,u)\right\|_{{\mathcal{H}}^{1}}&\lesssim\max_{j\in\\{1,\dots,m\\}}\left\|\hat{\chi}^{(3)}(\cdot,\omega)\right\|_{W^{1,\infty}(I_{j})}\left\||\tilde{u}|^{2}\tilde{u}+(\tilde{u}\cdot\tilde{u})\overline{\tilde{u}}\right\|_{{\mathcal{H}}^{1}}\\\
&\lesssim\|u\|^{3}_{{\mathcal{H}}^{1}}\end{split}$

Hence, $h(\cdot,\omega,u)\in{\mathcal{H}}^{1}$ for all $u\in
D(A)\cap{\mathcal{H}}^{1}$ and $\omega\in B_{\delta}(\omega_{0})$. For
$\varphi$, $\psi\in D(A)\cap{\mathcal{H}}^{1}$ we write
$\varphi=\left(\tilde{\varphi},\varphi_{3}\right)^{\mathsf{T}}$ with
$\tilde{\varphi}=(\varphi_{1},\varphi_{2})^{\mathsf{T}}$ and
$\psi=(\tilde{\psi},\psi_{3})^{\mathsf{T}}$ with
$\tilde{\psi}=\left(\psi_{1},\psi_{2}\right)^{\mathsf{T}}$ and obtain (94) as
follows:

$\begin{split}&\left\|h(\cdot,\omega,\psi)-h(\cdot,\omega,\varphi)\right\|_{{\mathcal{H}}^{1}}\\\
\lesssim&\max_{j\in\\{1,\dots,m\\}}\left\|\hat{\chi}^{(3)}(\cdot,\omega)\right\|_{W^{1,\infty}(I_{j})}\left(\left\||\tilde{\psi}|^{2}\tilde{\psi}-|\tilde{\varphi}|^{2}\tilde{\psi}\right\|_{{\mathcal{H}}^{1}}+\left\|(\tilde{\psi}\cdot\tilde{\psi})\overline{\tilde{\psi}}-(\tilde{\varphi}\cdot\tilde{\varphi})\overline{\tilde{\varphi}}\right\|_{{\mathcal{H}}^{1}}\right)\\\
\lesssim&\||\tilde{\psi}|^{2}\|_{{\mathcal{H}}^{1}}\left\|\tilde{\psi}-\tilde{\varphi}\right\|_{{\mathcal{H}}^{1}}+\left\||\tilde{\psi}|^{2}-|\tilde{\varphi}|^{2}\right\|_{{\mathcal{H}}^{1}}\left\|\tilde{\varphi}\right\|_{{\mathcal{H}}^{1}}+\left\|\tilde{\psi}\cdot\tilde{\psi}-\tilde{\varphi}\cdot\tilde{\varphi}\right\|_{{\mathcal{H}}^{1}}\left\|\tilde{\varphi}\right\|_{{\mathcal{H}}^{1}}\\\
\lesssim&\left(\left\|\psi\right\|^{2}_{{\mathcal{H}}^{1}}+\left\|\varphi\right\|^{2}_{{\mathcal{H}}^{1}}\right)\left\|\psi-\varphi\right\|_{{\mathcal{H}}^{1}}.\end{split}$

Finally, using the Lipschitz assumption on $\chi^{(3)}$, we get (95). Indeed,

$\begin{split}\left\|h(\cdot,\omega_{1},u)-h(\cdot,\omega_{2},u)\right\|_{{\mathcal{H}}^{1}}&\lesssim\max_{j\in\\{1,\dots,m\\}}\left\|\hat{\chi}^{(3)}_{1}(\cdot,\omega_{1})-\hat{\chi}^{(3)}_{2}(\cdot,\omega_{2})\right\|_{W^{1,\infty}(I_{j})}\left\||\tilde{u}|^{2}\tilde{u}+(\tilde{u}\cdot\tilde{u})\overline{\tilde{u}}\right\|_{{\mathcal{H}}^{1}}\\\
&\lesssim\left|\omega_{1}-\omega_{2}\right|\left\|u\right\|^{3}_{{\mathcal{H}}^{1}}.\end{split}$

∎

## Acknowledgment

Tomáš Dohnal acknowledges funding by the German Research Foundation (Deutsche
Forschungsgemeinschaft DFG), project-ID DO1467/4-1.

## References

  * [1] H. Alaeian and J. A. Dionne.  Parity-time-symmetric plasmonic metamaterials.  Phys. Rev. A, 89:033829, Mar 2014.
  * [2] J. N. Anker, W. P. Hall, O. Lyandres, N. C. Shah, J. Zhao, and R. P. Van Duyne.  Biosensing with plasmonic nanosensors.  Nature materials, 7(6):442–453, 2008.
  * [3] N. W. Ashcroft and N. D. Mermin.  Solid State Physics.  Holt-Saunders, 1976.
  * [4] W. L. Barnes, A. Dereux, and T. W. Ebbesen.  Surface plasmon subwavelength optics.  nature, 424(6950):824–830, 2003.
  * [5] D. Barton, M. Lawrence, H. Alaeian, B. Baum, and J. Dionne.  Parity-time symmetric plasmonics.  In D. Christodoulides and J. Yang, editors, Parity-time Symmetry and Its Applications, pages 301–349. Springer Singapore, Singapore, 2018.
  * [6] C. M. Bender and S. Boettcher.  Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry.  Phys. Rev. Lett., 80:5243–5246, Jun 1998.
  * [7] R. Boyd and D. Prato.  Nonlinear Optics.  Elsevier Science, 2008.
  * [8] M. Brown, T. Dohnal, M. Plum, and I. Wood.  Spectrum of the Maxwell Equations for a Flat Interface between Homogeneous Dispersive Media.  arXiv:2206.02037, 2022.
  * [9] A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar.  Self-focusing and spatial plasmon-polariton solitons.  Opt. Express, 17(24):21732–21737, 2009.
  * [10] T. Dohnal and G. Romani.  Eigenvalue bifurcation in doubly nonlinear problems with an application to surface plasmon polaritons.  Nonlinear Differential Equations and Applications NoDEA, 28(1):1–30, 2021.
  * [11] T. Dohnal and G. Romani.  Correction to: Eigenvalue bifurcation in doubly nonlinear problems with an application to surface plasmon polaritons.  Nonlinear Differential Equations and Applications NoDEA, 30(1):9, 2022.
  * [12] T. Dohnal and P. Siegl.  Bifurcation of eigenvalues in nonlinear problems with antilinear symmetry.  J. Math. Phys., 57(9):093502, 18, 2016.
  * [13] N. Dunford and J. Schwartz.  Linear Operators, Part 1: General Theory.  Wiley Classics Library. Wiley, 1988.
  * [14] E. Feigenbaum and M. Orenstein.  Plasmon-soliton.  Opt. Lett., 32(6):674–676, 2007.
  * [15] A. A. Govyadinov and V. A. Podolskiy.  Gain-assisted slow to superluminal group velocity manipulation in nanowaveguides.  Phys. Rev. Lett., 97:223902, Dec 2006.
  * [16] L. Guo, J. A. Jackman, H.-H. Yang, P. Chen, N.-J. Cho, and D.-H. Kim.  Strategies for enhancing the sensitivity of plasmonic nanosensors.  Nano Today, 10(2):213–239, 2015.
  * [17] J. Han, Y. Fan, L. Jin, Z. Zhang, Z. Wei, C. Wu, J. Qiu, H. Chen, Z. Wang, and H. Li.  Mode propagation in a pt-symmetric gain-metal-loss plasmonic system.  Journal of Optics, 16(4):045002, mar 2014.
  * [18] J. Homola, S. S. Yee, and G. Gauglitz.  Surface plasmon resonance sensors.  Sensors and actuators B: Chemical, 54(1-2):3–15, 1999.
  * [19] T. Kato.  Perturbation theory for linear operators, volume 132.  Springer Science & Business Media, 2013.
  * [20] S. Maier.  Plasmonics: Fundamentals and Applications.  Springer US, 2007.
  * [21] J. Mederski and W. Reichel.  Travelling waves for Maxwell’s equations in nonlinear and nonsymmetric media.  Nonlinear Differential Equations and Applications NoDEA, 30(2):22, 2022.
  * [22] H. Raether.  Surface Plasmons on Smooth and Rough Surfaces and on Gratings, volume 111.  Springer, 2006.
  * [23] R. H. Ritchie.  Plasma losses by fast electrons in thin films.  Phys. Rev., 106:874–881, Jun 1957.
  * [24] R. H. Ritchie.  Surface plasmons in solids.  Surface Science, 34(1):1–19, 1973.
  * [25] Y. Shen.  The Principles of Nonlinear Optics.  Pure & Applied Optics Series: 1-349. Wiley, 1984.
  * [26] I. I. Smolyaninov.  Quantum fluctuations of the refractive index near the interface between a metal and a nonlinear dielectric.  Phys. Rev. Lett., 94:057403, Feb 2005.
  * [27] C. A. Stuart.  Guidance properties of nonlinear planar waveguides.  Arch. Ration. Mech. Anal., 125(2):145–200, 1993.
  * [28] R. L. Sutherland, D. G. McLean, and S. Kirkpatrick.  Handbook of nonlinear optics.  Optical engineering. Marcel Dekker, New York, 2003.
  * [29] S. Tikhov and D. Valovik.  Theory of nonlinear transverse-magnetic guided waves in a plane waveguide filled with nonhomogeneous nonlinear medium.  Waves in Random and Complex Media, pages 1–25, 2022.
  * [30] V. K. Valev.  Characterization of nanostructured plasmonic surfaces with second harmonic generation.  Langmuir, 28(44):15454–15471, 2012.  PMID: 22889193.
  * [31] D. Valovik.  Novel propagation regimes for TE waves guided by a waveguide filled with Kerr medium.  Journal of Nonlinear Optical Physics & Materials, 25(04):1650051, 2016.
  * [32] D. V. Valovik and Y. G. Smirnov.  Propagation of TM waves in a Kerr nonlinear layer.  Computational Mathematics and Mathematical Physics, 48:2217–2225, 2008.