#  Bifurcations of thresholds in essential spectra of elliptic operators under
localized non-Hermitian perturbations

D. I. Borisov1, D. A. Zezyulin2111Corresponding author, M. Znojil3

###### Abstract

We consider the operator

$\mathcal{H}=\mathcal{H}^{\prime}-\frac{\partial^{2}\ }{\partial
x_{d}^{2}}\quad\text{on}\quad\omega\times\mathds{R}$

subject to the Dirichlet or Robin condition, where a domain
$\omega\subseteq\mathds{R}^{d-1}$ is bounded or unbounded. The symbol
$\mathcal{H}^{\prime}$ stands for a second order self-adjoint differential
operator on $\omega$ such that the spectrum of the operator
$\mathcal{H}^{\prime}$ contains several discrete eigenvalues $\Lambda_{j}$,
$j=1,\ldots,m$. These eigenvalues are thresholds in the essential spectrum of
the operator $\mathcal{H}$. We study how these thresholds bifurcate once we
add a small localized perturbation $\varepsilon\mathcal{L}(\varepsilon)$ to
the operator $\mathcal{H}$, where $\varepsilon$ is a small positive parameter
and $\mathcal{L}(\varepsilon)$ is an abstract, not necessarily symmetric
operator. We show that these thresholds bifurcate into eigenvalues and
resonances of the operator $\mathcal{H}$ in the vicinity of $\Lambda_{j}$ for
sufficiently small $\varepsilon$. We prove effective simple conditions
determining the existence of these resonances and eigenvalues and find the
leading terms of their asymptotic expansions. Our analysis applies to generic
non-self-adjoint perturbations and, in particular, to perturbations
characterized by the parity-time ($\mathcal{PT}$) symmetry. Potential
applications of our result embrace a broad class of physical systems governed
by dispersive or diffractive effects. As a case example, we employ our
findings to develop a scheme for a controllable generation of non-Hermitian
optical states with normalizable power and real part of the complex-valued
propagation constant lying in the continuum. The corresponding eigenfunctions
can be interpreted as an optical generalization of bound states embedded in
the continuum. For a particular example, the persistence of asymptotic
expansions is confirmed with direct numerical evaluation of the perturbed
spectrum.

1 Institute of Mathematics, Ufa Federal Research Center, Russian Academy of
Sciences, Ufa, Russia,
&
Bashkir State University, Ufa, Russia,
&
University of Hradec Králové, Hradec Králové, Czech Republic
<EMAIL_ADDRESS>
2 ITMO University, Saint-Petersburg 197101, Russia
<EMAIL_ADDRESS>
3 The Czech Academy of Sciences, Nuclear Physics Institute, Hlavní 130, 25068
Řež, Czech Republic
&
University of Hradec Králové, Hradec Králové, Czech Republic,
&
Institute of System Science, Durban University of Technology, P.O. Box 1334,
Durban 4000, South Africa
<EMAIL_ADDRESS>
##  1 Introduction

##### Physical context and motivation.

Physics of non-Hermitian Hamiltonians is attracting steadily growing attention
both on the fundamental level in the development of complex formulations of
quantum mechanics [2], [48], [47] and in several applied and experimental
fields, such as optics and photonics, Bose-Einstein condensates of atoms or
exciton-polaritons, acoustics, and in other areas where diffractive or
dispersive effects are governed by Schröginder-like elliptic operators, see
[70], [40], [24], [43], [22] for recent reviews. Prominent examples of
essentially non-Hermitian phenomena that were introduced in mathematical
literature long ago but entered various areas of physics much more recently
include exceptional points [35] and spectral singularities [53], [61], [64].
In particular, unusual effects associated with exceptional points are being
extensively discussed in optics and photonics (e.g. [30], [46], [57]), whereas
spectral singularities [53], [61], [64] are now understood to play an
important role in wave scattering [49], [50] and are used to implement
coherent perfect absorption of electromagnetic [66], sound [45], and matter
[52] waves. Another prototypical behavior, which is forbidden in Hermitian
physics but is of the utmost importance in non-Hermitian systems, is the
transition from the entirely real spectrum of eigenvalues to a partially
complex one.

In a real-world system, the non-Hermiticity usually corresponds to the
presence of an energy gain or absorption, which creates an effective complex
potential for propagating waves [51]. An especially interesting situation,
where a judicious balance between amplification and losses results in rich
physics, corresponds to the so-called parity-time ($\mathcal{PT}$) symmetric
systems famous for their property to robustly preserve reality of all
eigenvalues in spite of the absence of Hermiticity [3, 2]. Tuning a control
parameter of a $\mathcal{PT}$-symmetric system, one can realize various
qualitative changes in its spectral structure, and these changes are typically
associated with rich and intriguing behaviors. The simplest of those
behaviors, which was observed in a series of experiments [46, 57], corresponds
to the collision of a pair of real discrete eigenvalues in an exceptional
point with a subsequent splitting in a complex-conjugate pair. In systems with
a continuous spectrum, the situation is further enriched. In particular, the
bifurcation of an isolated eigenvalue from the bottom of the essential
spectrum can be accompanied by a so-called jamming anomaly, i.e., a non-
monotonous dependence of the energy flux through the gain-to-loss interface on
the parameter characterizing the strength of the non-Hermiticity [1]. The
bifurcations of a complex-conjugate pair of eigenvalues from an internal point
in the essential spectrum are even more interesting and have recently been
discussed as an unconventional mechanism of $\mathcal{PT}$-symmetry breaking
[67], [28], [41], [42] distinctively different from the better studied
$\mathcal{PT}$-symmetry breaking through an exceptional point. Complex
eigenfunctions associated with bifurcated eigenvalues are $L^{2}$-integrable,
and the real parts of these eigenvalues belong to the continuous spectrum.
This enables a physical interpretation of such eigenfunctions in terms of non-
Hermitian generalizations [44], [34] of bound states embedded in the
continuum, well-known in quantum mechanics [56], [63], [58], [60], optics, and
other fields [32]. It should be noticed at the same time that most of the
activity devoted to non-Hermitian optical bound states in the continuum is
being carried out for one-dimensional systems. For multi-dimensional
geometries, most of the available results are numerical in nature [34]. From
the practical point of view, it is also important that eigenfunctions
associated with emerging from the essential spectrum eigenvalues are extremely
weakly localized in the vicinity of the bifurcation, which hinders their
efficient numerical evaluation. Naturally, this problem is even more
pronounced in multi-dimensional geometries, where many more computational
resources are necessary to approximate the eigenfunctions. Therefore, any
analytical information on the properties of such states is highly desirable.

##### Mathematical context.

The phenomenon that a small localized perturbation of a self-adjoint
differential operator can generate discrete eigenvalues from the edges in the
essential spectrum is known for about a hundred of years. Its rigorous
mathematical study was initiated by classical works by B. Simon, M. Klaus, R.
Blankenbecler, M. L. Goldberger [62], [37], [4], [38] and since that time,
hundreds of papers on this subject were written. While classical works were
devoted to the Schrödinger operator on an axis and plane perturbed by a small
localized potentials, in further works the studies were made for plenty of
other models, like waveguide-like structures, see, for instance, [23], [21],
[54], [15], for periodic operators, see, for instance, [69], for operators
with distant perturbations [31] and many others. All these works treated
symmetric perturbations of self-adjoint operators, and the perturbed operators
were self-adjoint as well.

Non-symmetric perturbations of self-adjoint operators were studied in
essentially less details. In [25], [26] there was considered the Laplacian on
the axis perturbed by a small abstract localized operator, which was not
assumed to be symmetric. The main result was sufficient conditions ensuring
the existence and absence of the emerging eigenvalues from the bottom of the
essential spectrum and if they exist, the leading terms in the asymptotic
expansions for the emerging eigenvalues were found. These results were
essentially extended in [8], [9]. Here an unperturbed operator was an
arbitrary periodic self-adjoint operator on the line [8] or on the plane [9].
A perturbation was a small abstract operator not necessarily symmetric and
localized in a much weaker sense than in [25], [26]. The structure of the
spectra of such operators was studied in details. Qualitative properties like
stability of the essential spectrum, the countability of the point spectrum,
the absence of the residual spectrum, the existence of embedded eigenvalues
were addressed. Sufficient conditions ensuring the existence and absence of
the eigenvalues emerging from edges of internal gaps in the essential spectrum
were established and if they exist, the leading terms in their asymptotic
expansions were obtained. Eigenvalues emerging from the bottom of the
essential spectrum were also studied in [12], [13], [14] for waveguides with
$\mathcal{PT}$-symmetric Robin-type boundary condition. In [12], a planar
waveguide was considered with a locally perturbed coefficient in the
$\mathcal{PT}$-symmetric boundary condition. In [13], [14], similar two- and
three-dimensional waveguides were considered and the perturbation was a small
width of these waveguides. The main obtained results were sufficient
conditions ensuring the existence of the emerging eigenvalues and the leading
terms of their asymptotic expansions. We also mention work [27], where the
Dirichlet or Neumann Laplacian in a multi-dimensional cylinder was considered
and it was perturbed by a small localized non-symmetric perturbations. The
eigenvalues bifurcating from the bottom and the internal thresholds in the
essential spectrum were studied. There were obtained certain sufficient
conditions ensuring the existence of such eigenvalues and the leading terms of
their asymptotic expansions were calculated. However, there was a gap in
calculations in [27], which made the results of this work true but incomplete.
Namely, while working with operators providing meromorphic continuations for
the resolvent in the vicinity of internal thresholds, the author of [27]
considered only one continuation, while, as we show in the present work, even
in our more general setting two continuations exist and complex eigenvalues
are the poles just for one of these continuations. This is why the results of
[27] described, roughly speaking, only half of emerging eigenvalues.

The bifurcations of the thresholds in the essential spectrum can be also
studied for perturbations of non-self-adjoint operators provided the spectral
structure of the limiting operator is known in sufficiently great details. As
examples, we mention works [20], [65], where an evolutionary nonlinear
Schrödinger equation was considered with both linear and nonlinear
perturbations. The linearization of this equation on solitary wave solutions
gave rise to a spectral problem for a linear non-self-adjoint operator. The
essential spectral of such operator consists of two real semi-axes; the main
results of [20], [65] provided conditions, under which the end-points of the
essential spectrum bifurcated into eigenvalues. If the latter existed, their
two-terms asymptotic expansions were calculated.

An important feature of the eigenvalues emergence is that usually the total
multiplicity of the emerging eigenvalues does not exceed the multiplicity of
edge in the essential spectrum from which they emerge. The multiplicity of the
edge is to be treated in the sense of some appropriate generalized
eigenfunctions. However, there were found examples, when this commonly
believed rule failed. The earliest work on this subject we know is paper [29],
where the Schrödinger operator on $\mathds{R}^{3}$ perturbed by a small
localized potential was considered. It was found that in certain cases, an
$n$-multiple bottom of the essential spectrum can generate $n$ eigenvalues and
$n$ anti-bound states or $2n$ resonances. In [16], a similar phenomenon was
found for the Dirichlet Laplacian in a pair of three-dimensional layers
coupled by a window, when the perturbation was a small variation of the window
shape. Very recently we succeeded to find an even more impressive example of
infinitely many eigenvalues and/or resonances emerging from the bottom of an
essential spectrum. This was done in papers [17], [18], where we considered a
one-dimensional Schrödinger operator on the axis with two complex localized
potentials, the supports of which were separated by a large distance. It was
found that as this distance increases, more and more resonances and
eigenvalues appear in the vicinity of the bottom of the essential spectrum,
while the multiplicity of this bottom is at most one. The location and
asymptotic behaviour of these emerging eigenvalues and resonances were
analyzed in details.

Emerging eigenvalues were also studied not only in the context of classical
eigenvalue problems, but also for more complicated operator pencils. In [11],
there was considered a special quadratic operator pencil on the line with a
special small periodic $\mathcal{PT}$-symmetric perturbation. The structure of
the gaps in the essential spectrum and complex eigenvalues in the vicinities
of the edges of these gaps were analyzed in great details. In [10], there was
considered a similar quadratic operator pencil with a special small localized
$\mathcal{PT}$-symmetric potential and there were studied eigenvalues emerging
from thresholds in the essential spectrum. Sufficient existence conditions
were established and the leading terms of the asymptotic expansions of the
emerging eigenvalues were obtained.

##### Model and results.

In the present paper we carry out a rigorous analysis of bifurcations of
isolated eigenvalues and resonances from the essential spectrum of a multi-
dimensional operator under a small localized general abstract perturbation.
Namely, we consider a self-adjoint operator of the form

$\mathcal{H}=\mathcal{H}^{\prime}-\frac{\partial^{2}\ }{\partial
x_{d}^{2}}\quad\text{on}\quad\omega\times\mathds{R}$

subject to the Dirichlet or Robin condition, where
$\omega\subseteq\mathds{R}^{d-1}$ is some domain, which can be both bounded
and unbounded and $\mathcal{H}^{\prime}$ is a self-adjoint second order
differential operator on $\omega$ subject to the same boundary condition as
$\mathcal{H}$. We assume that the spectrum of the operator
$\mathcal{H}^{\prime}$ contains several discrete eigenvalues
$\Lambda_{1}\leqslant\Lambda_{2}\leqslant\ldots\leqslant\Lambda_{m}$ below the
essential spectrum. Then the essential spectrum of the operator $\mathcal{H}$
is the half-line $[\Lambda_{1},+\infty)$ and the mentioned eigenvalues become
thresholds in this essential spectrum. We add a small localized perturbation
to the operator $\mathcal{H}$. This perturbation reads as
$\varepsilon\mathcal{L}(\varepsilon)$, where $\varepsilon$ is a small positive
parameter and $\mathcal{L}(\varepsilon)$ is an abstract not necessarily
symmetric operator acting from a weighted Sobolev space
$W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ into a weighted Lebesgue space
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$. Exact definitions of these spaces
will be given in the next section and now we just say that these weights and
the operator $\mathcal{L}(\varepsilon)$ are designed so that this operator
maps exponentially growing functions into exponentially decaying ones. The
latter fact is exactly how we understand the localization of this operator in
a generalized sense.

The main result of our paper describes how the thresholds $\Lambda_{j}$ in the
essential spectrum of the operator $\mathcal{H}$ bifurcate under the presence
of the perturbation $\varepsilon\mathcal{L}(\varepsilon)$. We show that if the
bottom of the essential spectrum $\Lambda_{1}$ is an $m$-multiple eigenvalue
of the operator $\mathcal{H}^{\prime}$, then there can be at most $m$
eigenvalues and resonances of the operator $\mathcal{H}$ in the vicinity of
$\Lambda_{1}$ for sufficiently small $\varepsilon$. The vicinity of an
internal threshold $\Lambda_{j}>\Lambda_{1}$ in the essential spectrum being
an $m$-multiple eigenvalue of the operator $\mathcal{H}^{\prime}$ can contain
at most $2m$ eigenvalues and resonances of the operator $\mathcal{H}$. The
eigenvalues and resonances are identified via an analysis of the poles of an
appropriate meromorphic continuation of the resolvent in the vicinity of each
threshold $\Lambda_{j}$, $j\geqslant 1$. Each such pole generates either an
eigenvalue or a resonance, and we provide simple sufficient conditions
allowing one to identify whether a considered pole is an eigenvalue or a
resonance. We also construct two-terms asymptotic expansions for the emerging
eigenvalues and resonances.

##### Applications to specific models.

While our result is rather general and applies to a broad range of physical
models, where elliptic operators play the prominent role, we will exemplify
applications of the work using some particular physical models. Namely, we
discuss various examples of the unperturbed operator $\mathcal{H}$ and of the
perturbation $\mathcal{L}(\varepsilon)$. Then we consider models of two- and
three-dimensional waveguides and models of two- and three-dimensional quantum
oscillators. As a perturbation, we choose a small localized complex potential.
Such choice is motivated by physical models of optical waveguides filled with
a homogeneous medium, when the refractive index of the waveguide is locally
modulated. This creates a small, generally, non-Hermitian perturbation in the
form of an effective complex-valued optical potential. In particular, this
potential can be $\mathcal{PT}$-symmetric. Another physical model motivating
the above examples is a two-dimensional Bose-Einstein condensate trapped in a
harmonic potential in one dimension and without any trapping in the second
dimension. The nonlinear interactions between particles of the condensate are
assumed to be negligible such that its evolution can be described by the
linear Schrödinger operator. As a perturbation, a localized non-Hermitian
defect such as a localized dissipation serves. The similar approach can be
applied to a three-dimensional condensate, where a localized perturbation can
trigger formation of fully localized structures with internal vorticity.

For such examples we show that given an internal threshold of a multiplicity
$n$, by tuning appropriately the perturbing potential, we can make the
threshold to bifurcate into $n$ pairs of complex-conjugated eigenvalues. This
example demonstrates that the total multiplicity of the emerging eigenvalues
can exceed the multiplicity of the internal threshold.

##### Organization of the paper.

The rest of this paper is organized as follows. In Section 2 we elaborate
rigorous mathematical formulation of the problem, and then present and discuss
the main results which are formulated in several theorems. Section 3 is
dedicated to examples, including a case study of optical bound states in the
continuum emerging under a small $\mathcal{PT}$-symmetric perturbation.
Sections 4 and 5 contain the proofs of theorems.

##  2 Problem and results

###  2.1 Problem

Let $x^{\prime}=(x_{1},\ldots,x_{d-1})$, $x=(x^{\prime},x_{d})$ be Cartesian
coordinates in $\mathds{R}^{d-1}$ and $\mathds{R}^{d}$, respectively, where
$d\geqslant 2$, and $\omega\subseteq\mathds{R}^{d-1}$ be an arbitrary domain.
The domain $\omega$ can be bounded or unbounded, the case
$\omega=\mathds{R}^{d-1}$ is also possible. If the boundary of the domain
$\omega$ is non-empty, we assume that $\partial\omega\in C^{2}$. We let
$\Omega:=\omega\times\mathds{R}$ and we suppose that the domain $\omega$ is
such that

$\|u\|_{L_{2}(\partial\omega)}\leqslant C\|u\|_{W_{2}^{1}(\omega)}$

for all $u\in W_{2}^{1}(\omega)$ with a constant $C$ independent of $u$. This
inequality implies that

$\|u\|_{L_{2}(\partial\Omega)}\leqslant C\|u\|_{W_{2}^{1}(\Omega)}$

for all $u\in W_{2}^{1}(\Omega)$ with a constant $C$ independent of $u$. This
means that on the boundary of the domain $\Omega$, the traces of the functions
in $W_{2}^{1}(\Omega)$ are well-defined and the trace operator is bounded.
This fact is employed below in definitions of various operators and
sesquilinear forms without explicit mentioning.

By $A_{ij}=A_{ij}(x^{\prime})$, $A_{j}=A_{j}(x^{\prime})$,
$A_{0}=A_{0}(x^{\prime})$, $i,j=1,\ldots,d-1$, we denote real functions
defined on $\overline{\omega}$ and with the following smoothness:
$A_{ij},A_{j}\in C^{1}(\overline{\omega})$, $A_{0}\in C(\overline{\omega})$.
The functions $A_{ij}$ satisfy the usual uniform ellipticity condition, that
is, $A_{ij}=A_{ji}$ and

$\sum\limits_{i,j=1}^{d-1}A_{ij}{(x^{\prime})}\xi_{i}\overline{\xi_{j}}\geqslant
c_{0}\sum\limits_{i=1}^{d-1}|\xi_{i}|^{2}\quad\text{for all}\quad
x^{\prime}\in\overline{\omega},\quad\xi_{i}\in\mathds{C},$

where $c_{0}>0$ is a positive constant independent of $x^{\prime}$ and
$\xi_{i}$. The functions $A_{ij}$ and $A_{j}$ are assumed to be uniformly
bounded on $\overline{\omega}$, while for $A_{0}$ only an uniform lower bound
is supposed. By $\mathrm{i}$ we denote the imaginary unit.

In terms of the introduced functions we define an operator

$\mathcal{H}=-\sum\limits_{i,j=1}^{d-1}\frac{\partial\ }{\partial
x_{i}}A_{ij}\frac{\partial\ }{\partial x_{j}}-\frac{\partial^{2}\ }{\partial
x_{d}^{2}}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial\
}{\partial x_{j}}+\frac{\partial\ }{\partial
x_{j}}A_{j}\right)+A_{0}\quad\text{in}\quad\Omega$ (2.1)

subject to the Dirichlet condition or Robin condition:

$u=0\quad\text{on}\quad\partial\Omega\qquad\text{or}\qquad\frac{\partial
u}{\partial\boldsymbol{\nu}}-au=0\quad\text{on}\quad\partial\Omega.$ (2.2)

In the case of Robin condition, the conormal derivative is defined as

$\frac{\partial
u}{\partial\boldsymbol{\nu}}:=\sum\limits_{i,j=1}^{d-1}A_{ij}\nu_{i}\frac{\partial
u}{\partial
x_{j}}-\mathrm{i}\sum\limits_{j=1}^{d-1}A_{j}\nu_{j}{u}+\nu_{d}\frac{\partial
u}{\partial x_{d}},$

where $\nu=(\nu_{1},\ldots,\nu_{d})$ is the unit outward normal to
$\partial\Omega$ and $a=a(x^{\prime})$ is a real function defined on
$\partial\Omega$. We assume that $a\in C(\partial\Omega)$ and that this
function is uniformly bounded on $\partial\Omega$. We define a chosen boundary
operator in (2.2) by $\mathcal{B}$, that is, $\mathcal{B}u=u$ or
$\mathcal{B}u=\frac{\partial u}{\partial\boldsymbol{\nu}}-au$.

Rigorously we introduce the operator $\mathcal{H}$ as follows. In the space
$L_{2}(\Omega)$ we define a sesquilinear form

$\displaystyle\mathfrak{h}(u,v):=$
$\displaystyle\sum\limits_{i,j=1}^{d-1}\left(A_{ij}\frac{\partial u}{\partial
x_{j}},\frac{\partial v}{\partial
x_{i}}\right)_{L_{2}(\Omega)}+\left(\frac{\partial u}{\partial
x_{d}},\frac{\partial v}{\partial
x_{d}}\right)_{L_{2}(\Omega)}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial
u}{\partial x_{j}},v\right)_{L_{2}(\Omega)}$
$\displaystyle-\mathrm{i}\sum\limits_{j=1}^{d-1}\left(u,A_{j}\frac{\partial
v}{\partial x_{j}}\right)_{L_{2}(\Omega)}+(A_{0}u,v)_{L_{2}(\Omega)}$

on the domain $\mathfrak{D}(\mathfrak{h}):=\mathring{W}_{2}^{1}(\Omega)\cap
L_{2}(\Omega,(1+|A_{0}|)dx)$ if the Dirichlet condition is chosen in (2.2),
and

$\displaystyle\mathfrak{h}(u,v):=$
$\displaystyle\sum\limits_{i,j=1}^{d-1}\left(A_{ij}\frac{\partial u}{\partial
x_{j}},\frac{\partial v}{\partial
x_{i}}\right)_{L_{2}(\Omega)}+\left(\frac{\partial u}{\partial
x_{d}},\frac{\partial v}{\partial
x_{d}}\right)_{L_{2}(\Omega)}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial
u}{\partial x_{j}},v\right)_{L_{2}(\Omega)}$
$\displaystyle-\mathrm{i}\sum\limits_{j=1}^{d-1}\left(u,A_{j}\frac{\partial
v}{\partial
x_{j}}\right)_{L_{2}(\Omega)}+(A_{0}u,v)_{L_{2}(\Omega)}-(au,v)_{L_{2}(\partial\Omega)}$

on the domain $\mathfrak{D}(\mathfrak{h}):=W_{2}^{1}(\Omega)\cap
L_{2}(\Omega,(1+|A_{0}|)dx)$. Here $\mathring{W}_{2}^{1}(\Omega)$ is a
subspace of the space $W_{2}^{1}(\Omega)$ consisting of the functions with a
zero trace on $\partial\Omega$. Given a positive function $\phi$ on $\Omega$,
by $L_{2}(\Omega,\phi dx)$ we denote a weighted space formed by the functions
in $L_{2,loc}(\Omega)$ with a finite norm $\|\cdot\|_{L_{2}(\Omega,\phi dx)}$
defined as

$\|u\|_{L_{2}(\Omega,\phi\,dx)}^{2}=\int\limits_{\Omega}|u(x)|^{2}\phi\,dx.$

Thanks to the above assumptions on the functions $A_{ij}$, $A_{j}$, $A_{0}$
and $a$, the form $\mathfrak{h}$ is closed, symmetric and lower-semibounded.
The self-adjoint operator in $L_{2}(\Omega)$ associated with this form is
exactly the operator $\mathcal{H}$.

We introduce one more weighted space
$W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ as a subspace of
$W_{2,loc}^{2}(\Omega)$ formed by the functions with finite norms
$\|\cdot\|_{W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)}$, is defined as
follows:

$\|u\|_{W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)}^{2}=\int\limits_{\Omega}\sum\limits_{\begin{subarray}{c}\alpha\in\mathds{Z}_{+}^{2}\\\
|\alpha|\leqslant
2\end{subarray}}|\partial^{\alpha}u(x)|^{2}e^{-{\vartheta}|x_{d}|}dx.$

Here $\vartheta>0$ is some fixed constant. By $\varepsilon$ we denote a small
positive parameter and the symbols $\mathcal{L}_{1}$, $\mathcal{L}_{2}$,
$\mathcal{L}_{3}=\mathcal{L}_{3}(\varepsilon)$ stand for operators mapping the
space $W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ into
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$. These operators are assumed to be
bounded; the operator $\mathcal{L}_{3}$ is bounded uniformly in $\varepsilon$.
We stress that the operators $\mathcal{L}_{1}$, $\mathcal{L}_{2}$,
$\mathcal{L}_{3}$ are not supposed to be necessarily symmetric.

The main object of our study is a perturbed operator

$\mathcal{H}_{\varepsilon}=\mathcal{H}+\varepsilon\mathcal{L}(\varepsilon),\qquad\mathcal{L}(\varepsilon):=\mathcal{L}_{1}+\varepsilon\mathcal{L}_{2}+\varepsilon^{2}\mathcal{L}_{3}(\varepsilon)$

on $\mathfrak{D}(\mathcal{H})$. The operator $\mathcal{H}_{\varepsilon}$ is
well-defined since

$\mathfrak{D}(\mathcal{H})\subseteq W_{2}^{2}(\Omega)\subseteq
W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx).$

Moreover, it is clear that the operator $\mathcal{L}(\varepsilon)$ is
relatively bounded with respect to the operator $\mathcal{H}$ and this is why,
for sufficiently small $\varepsilon$, the operator $\mathcal{H}_{\varepsilon}$
is closed.

Our main aim is to study the behaviour of the eigenvalues of the operator
$\mathcal{H}_{\varepsilon}$ emerging from certain internal points in its
essential spectrum. We denote the latter by
$\operatorname{\sigma_{ess}}(\cdot)$ and define it in terms of a
characteristic sequences. Namely, a point $\lambda$ belongs to an essential
spectrum $\operatorname{\sigma_{ess}}(\mathcal{A})$ of some operator
$\mathcal{A}$ if there exists a bounded noncompact sequence
$u_{d}\in\mathfrak{D}(\mathcal{A})$ such that

$\inf\limits_{d}\|u_{d}\|>0\qquad\text{and}\qquad(\mathcal{A}-\lambda)u_{d}\to
0,\quad d\to\infty.$

In order to describe the essential spectrum
$\operatorname{\sigma_{ess}}(\mathcal{H}_{\varepsilon})$, we introduce two
auxiliary operators $\mathcal{H}^{\prime}$ and $\mathcal{H}_{0}$. The former
is a self-adjoint operator in $L_{2}(\mathds{R}^{d-1})$ associated with a
lower-semibounded symmetric sesquilinear form

$\displaystyle\mathfrak{h}^{\prime}(u,v):=$
$\displaystyle\sum\limits_{i,j=1}^{d-1}\left(A_{ij}\frac{\partial u}{\partial
x_{j}},\frac{\partial v}{\partial
x_{i}}\right)_{L_{2}(\omega)}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial
u}{\partial x_{j}},v\right)_{L_{2}(\omega)}$
$\displaystyle-\mathrm{i}\sum\limits_{j=1}^{d-1}\left(u,A_{j}\frac{\partial
v}{\partial x_{j}}\right)_{L_{2}(\omega)}+(A_{0}u,v)_{L_{2}(\omega)}$

on the domain
$\mathfrak{D}(\mathfrak{h}^{\prime}):=\mathring{W}_{2}^{1}(\omega)\cap
L_{2}(\omega,(1+|A_{0}|)dx^{\prime})$ if the Dirichlet condition is chosen in
(2.2) and

$\displaystyle\mathfrak{h}^{\prime}(u,v):=$
$\displaystyle\sum\limits_{i,j=1}^{d-1}\left(A_{ij}\frac{\partial u}{\partial
x_{j}},\frac{\partial v}{\partial
x_{i}}\right)_{L_{2}(\omega)}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial
u}{\partial x_{j}},v\right)_{L_{2}(\omega)}$
$\displaystyle-\mathrm{i}\sum\limits_{j=1}^{d-1}\left(u,A_{j}\frac{\partial
v}{\partial
x_{j}}\right)_{L_{2}(\omega)}+(A_{0}u,v)_{L_{2}(\omega)}-(au,v)_{L_{2}(\partial\omega)}$

on the domain $\mathfrak{D}(\mathfrak{h}^{\prime}):=W_{2}^{1}(\omega)\cap
L_{2}(\omega,(1+|A_{0}|)dx^{\prime})$ if the Robin condition is chosen in
(2.2). This is the operator

$\mathcal{H}^{\prime}=-\sum\limits_{i,j=1}^{d-1}\frac{\partial\ }{\partial
x_{i}}A_{ij}\frac{\partial\ }{\partial
x_{j}}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial\ }{\partial
x_{j}}+\frac{\partial\ }{\partial
x_{j}}A_{j}\right)+A_{0}\quad\text{in}\quad\omega$

subject to the Dirichlet condition or Robin condition:

$u=0\quad\text{on}\quad\partial\omega\qquad\text{or}\qquad\frac{\partial
u}{\partial\boldsymbol{\nu}^{\prime}}=au\quad\text{on}\quad\partial\Omega,\qquad\frac{\partial
u}{\partial\boldsymbol{\nu}^{\prime}}:=\sum\limits_{i,j=1}^{d-1}A_{ij}\nu_{i}\frac{\partial
u}{\partial x_{j}}-\mathrm{i}\sum\limits_{j=1}^{d-1}A_{j}\nu_{j}.$

The operator $\mathcal{H}_{0}$ is a one-dimensional Schrödinger operator

$\mathcal{H}_{0}:=-\frac{d^{2}\ }{dx_{d}^{2}}$

in $L_{2}(\mathds{R})$ on the domain $W_{2}^{2}(\mathds{R})$. Its spectrum is
pure essential and coincides with $[0,+\infty)$. We assume that there exists a
constant $c_{0}$ such that the spectrum of the operator $\mathcal{H}^{\prime}$
below this constant consists of finitely many discrete eigenvalues, which we
denote by $\Lambda_{j}$ and we arrange them in an ascending order counting
multiplicities:

$\Lambda_{1}\leqslant\Lambda_{2}\leqslant\ldots\leqslant\Lambda_{m}<c_{0}.$

The associated orthonormalized in $L_{2}(\mathds{R}^{d-1})$ eigenfunctions are
denoted by $\psi_{j}=\psi_{j}(x^{\prime})$, $j=1,\ldots,m$.

The essential spectrum of the operator $\mathcal{H}_{\varepsilon}$ is
described in the following lemma.

######  Lemma 2.1.

The essential spectrum of the operator $\mathcal{H}_{\varepsilon}$ coincides
with that of the operator $\mathcal{H}$ for all sufficiently small
$\varepsilon$ and is given by the identity:

$\operatorname{\sigma_{ess}}(\mathcal{H}_{\varepsilon})=\operatorname{\sigma_{ess}}(\mathcal{H})=\operatorname{\sigma}(\mathcal{H})=[\Lambda_{1},+\infty).$

where $\operatorname{\sigma}(\cdot)$ denotes a spectrum of an operator.

According to this lemma, the points $\Lambda_{j}$, $j=1,\ldots,m$, belong to
the essential spectrum of the operator $\mathcal{H}_{\varepsilon}$. The point
$\Lambda_{1}$ is the bottom of such spectrum, while other points $\Lambda_{j}$
are internal thresholds.

###  2.2 Main results

Our results describe a meromorphic continuation of the resolvent of the
operator $\mathcal{H}_{\varepsilon}$ in the vicinity of the points
$\Lambda_{j}$ as well as eigenvalues and resonances emerging from these points
due to the presence of the perturbation $\varepsilon\mathcal{L}(\varepsilon)$.
Before presenting our main results, we introduce some auxiliary constants and
notations.

By $B_{\delta}$ we denote a ball of radius $\delta$ centered at the origin in
the complex plane. We fix $p\in\\{1,\ldots,m\\}$ and assume that
$\Lambda_{p}=\ldots=\Lambda_{p+n-1}$ is an $n$-multiple eigenvalue of the
operator $\mathcal{H}^{\prime}$, where $n\geqslant 1$. Then we consider a new
complex parameter $k$ ranging in a small neighbourhood of the origin and we
introduce auxiliary functions:

$\displaystyle
K_{j}(k):=-\mathrm{i}\sqrt{\Lambda_{p}-\Lambda_{j}-k^{2}}\qquad\text{as}\quad
j<p,$ $\displaystyle
K_{j}(k):=k\hphantom{1\sqrt{\Lambda_{p}-\Lambda_{j}-k^{2}}}\qquad\text{as}\quad
j=p,\ldots,p+n-1,$ $\displaystyle
K_{j}(k):=\sqrt{\Lambda_{j}-\Lambda_{p}+k^{2}}\hphantom{-\mathrm{i}}\qquad\text{as}\quad
j\geqslant p+n.$

Hereinafter the branch of the square root is fixed by the condition
$\sqrt{1}=1$ with the branch cut along the negative real semi-axis. Given
$R>0$, we let $\Omega_{\pm}^{R}:=\Omega\cap\\{x:\pm x_{d}>R\\}$.

Now we are in position to formulate our first main result.

######  Theorem 2.1.

Fix $p\in\\{1,\ldots,m\\}$, and $\tau\in\\{-1,+1\\}$ and let
$\Lambda_{p}=\ldots=\Lambda_{p+n-1}$ be an $n$-multiple eigenvalue of the
operator $\mathcal{H}^{\prime}$, where $n\geqslant 1$. For all sufficiently
small $\varepsilon$, the resolvent
$(\mathcal{H}_{\varepsilon}{-}\Lambda_{p}+k^{2})^{-1}$ admits a meromorphic
continuation with respect to a complex parameter $k$ ranging in a sufficiently
small neighbourhood of the origin. Namely, there exists a bounded operator

$\mathcal{R}_{\varepsilon,\tau}(k):\,L_{2}(\Omega,e^{{\vartheta}|x_{d}|}\,dx)\to
W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}\,dx)$

meromorphic with respect to complex $k\in B_{\delta}$ for a sufficiently small
fixed $\delta$ independent of $\varepsilon$. If $p=1$, then the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ is independent of the choice of $\tau$ and
for $\operatorname{Re}k>0$, this operator coincides with the resolvent
$(\mathcal{H}_{\varepsilon}-\Lambda_{1}+k^{2})^{-1}$ restricted on
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}\,dx)$. If $p>1$, then the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ does depend on the choice of $\tau$ and
for $\operatorname{Re}k>0$ and $\tau\operatorname{Im}k^{2}<0$, this operator
coincides with the resolvent
$(\mathcal{H}_{\varepsilon}-\Lambda_{p}+k^{2})^{-1}$ restricted on
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}\,dx)$.

For all $f\in L_{2}(\Omega,e^{{\vartheta}|x_{d}|}\,dx)$, the function
$u_{\varepsilon}:=\mathcal{R}_{\varepsilon,\tau}(k)f$ solves the boundary
value problem

$\begin{gathered}\left(-\sum\limits_{i,j=1}^{d-1}\frac{\partial\ }{\partial
x_{i}}A_{ij}\frac{\partial\ }{\partial x_{j}}-\frac{\partial^{2}\ }{\partial
x_{d}^{2}}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial\
}{\partial x_{j}}+\frac{\partial\ }{\partial
x_{j}}A_{j}\right)+A_{0}+\varepsilon\mathcal{L}(\varepsilon)-\Lambda_{p}+k^{2}\right)u_{\varepsilon}=f\quad\text{in}\quad\Omega,\\\
\mathcal{B}u_{\varepsilon}=0\quad\text{on}\quad\partial\Omega,\end{gathered}$
(2.3)

and for sufficiently large $x_{d}$ it can be represented as follows:

$u_{\varepsilon}(x,k)=\sum\limits_{j=1}^{m}u_{\varepsilon,j}^{\pm}(x_{d},k)\psi_{j}(x^{\prime})+u_{\varepsilon,\bot}^{\pm}(x,k),\qquad\pm
x_{d}>R,$ (2.4)

where $R$ is some fixed number, $u_{\varepsilon,j}^{\pm}\in
L_{2}(I_{\pm},e^{{\mp\vartheta}x_{d}}dx_{d})$ are some meromorphic in $k\in
B_{\delta}$ functions, $I_{+}:=(R,+\infty)$, $I_{-}:=(-\infty,-R)$, possessing
the asymptotic behavior

$\displaystyle u_{\varepsilon,j}^{\pm}(x_{d})=e^{-\tau
K_{j}(k)|x_{d}|}\big{(}C_{\varepsilon,j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},\quad$
$\displaystyle|x_{d}|\to\infty,\qquad j=1,\ldots,p-n+1,$ (2.5) $\displaystyle
u_{\varepsilon,j}^{\pm}(x_{d})=e^{-K_{j}(k)|x_{d}|}\big{(}C_{\varepsilon,j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},$
$\displaystyle|x_{d}|\to\infty,\qquad j=p,\ldots,m,$

$C_{\varepsilon,j}^{\pm}(k)$ are some meromorphic in $k\in B_{\delta}$
functions, $0<{\tilde{\vartheta}<\vartheta}$ is some fixed constant
independent of $k$ and $x$, and $u_{\varepsilon,{\bot}}^{\pm}\in
W_{2}^{2}(\Omega_{R}^{\pm})$ are some functions meromorphic in $k\in
B_{\delta}$ and obeying the identities

$(u_{\varepsilon,\bot}^{\pm}(\cdot,x_{d}),\psi_{j})_{L_{2}({\omega})}=0$ (2.6)

for almost each $x_{d}\in I_{\pm}$ and for each $j=1,\ldots,m$.

If $k_{\varepsilon}\in B_{\delta}$ is a pole of the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$, for $k=k_{\varepsilon}$, problem (2.3)
with $f=0$ has a non-trivial solution $\psi_{\varepsilon}$ in
$W_{2,loc}^{2}(\Omega)$, which satisfies a representation similar to (2.4):

$\psi_{\varepsilon}(x)=\sum\limits_{j=1}^{m}\phi_{\varepsilon,j}^{\pm}(x_{d},k)\psi_{j}(x^{\prime})+\psi_{\varepsilon,\bot}^{\pm}(x,k),\qquad\pm
x_{d}>R,$ (2.7)

where $R$ is some fixed number, $\phi_{\varepsilon,j}^{\pm}\in
L_{2}(I_{\pm},e^{{\mp\vartheta}x_{d}}dx_{d})$ are functions with the
asymptotic behavior

$\displaystyle\phi_{\varepsilon,j}^{\pm}(x_{d})=e^{-\tau
K_{j}(k)|x_{d}|}\big{(}c_{\varepsilon,j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},\quad$
$\displaystyle|x_{d}|\to\infty,\qquad j=1,\ldots,p-n+1,$ (2.8)
$\displaystyle\phi_{\varepsilon,j}^{\pm}(x_{d})=e^{-K_{j}(k)|x_{d}|}\big{(}c_{\varepsilon,j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},$
$\displaystyle|x_{d}|\to\infty,\qquad j=p,\ldots,m,$

$c_{\varepsilon,j}^{\pm}(k)$ are some constants, and
$\psi_{\varepsilon}^{\pm}\in W_{2}^{2}(\Omega_{R}^{\pm})$ are some functions
obeying the identities

$(\psi_{\varepsilon,\bot}^{\pm}(\cdot,x_{d}),\psi_{j})_{L_{2}({\omega})}=0$
(2.9)

for almost each $x_{d}\in I_{\pm}$ and for each $j=1,\ldots,m$.

We define a subspace $L^{\bot}$ in $L_{2}(\Omega)$ as a set of functions $v\in
L_{2}(\Omega)$ such that

$(v(\cdot,x_{d}),\psi_{j})_{L_{2}(\omega)}=0$

for almost each $x_{d}\in\mathds{R}$ and for all $j=1,\ldots,m$. The space
$L^{\bot}$ is a Hilbert one. By $\mathcal{H}^{\bot}$ we denote the restriction
of the operator $\mathcal{H}$ on $\mathfrak{D}(\mathcal{H})\cap L^{\bot}$. The
following lemma will be proved in Section 4.1.

######  Lemma 2.2.

The space $L^{\bot}$ is invariant for the operator $\mathcal{H}^{\bot}$, that
is, this operator maps $\mathfrak{D}(\mathcal{H})\cap L^{\bot}$ into
$L^{\bot}$. This is an unbounded self-adjoint operator in $L^{\bot}$ and its
spectrum is located in $[c_{0},+\infty)$.

The above lemma means that the resolvent
$(\mathcal{H}^{\bot}-\Lambda_{p})^{-1}$ is well-defined for all $p=1,\ldots,m$
as an operator from $L_{2}(L^{\bot})$ into $\mathfrak{D}(\mathcal{H})\cap
L^{\bot}$. As above, we fix $p\in\\{1,\ldots,m\\}$ and assume that
$\Lambda_{p}=\ldots=\Lambda_{p+n-1}$, where $n\geqslant 1$, and in terms of
the latter resolvent, we introduce an auxiliary operator mapping
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ into
$W_{2}^{2}(\Omega,e^{{-\vartheta}|x_{d}|}dx)$:

$\displaystyle\begin{aligned}
(\mathcal{G}_{p,\tau}f)(x):=&\sum\limits_{j=1}^{p-1}\frac{\psi_{j}(x^{\prime})}{2\tau
K_{j}(0)}\int\limits_{\Omega}e^{-\tau
K_{j}(0)|x_{d}-y_{d}|}\overline{\psi_{j}(y^{\prime})}f(y)\,dy\\\
&-\frac{1}{2}\sum\limits_{j=p}^{p+n-1}\psi_{j}(x^{\prime})\int\limits_{\Omega}|x_{d}-y_{d}|\overline{\psi_{j}(y^{\prime})}f(y)\,dy\\\
&+\sum\limits_{j=p+n}^{m}\frac{\psi_{j}(x^{\prime})}{2K_{j}(0)}\int\limits_{\Omega}e^{-K_{j}(0)|x_{d}-y_{d}|}\overline{\psi_{j}(y^{\prime})}f(y)\,dy+((\mathcal{H}^{\bot}-\Lambda_{p})^{-1}f^{\bot})(x),\end{aligned}$
(2.10) $\displaystyle
f^{\bot}(x):=f(x)-\sum\limits_{j=1}^{m}f_{j}(x_{d})\psi_{j}(x^{\prime}).$
(2.11)

As above, here $\tau\in\\{-1,+1\\}$. In the case $p=1$, the first sum in the
above definition is missing and the operator $\mathcal{G}_{p,\tau}$ becomes
independent of the choice of $\tau$.

We define the matrix $\mathrm{M}_{1}$ with entries

$M_{1}^{ij}:=-\frac{1}{2}\int\limits_{\Omega}\overline{\psi_{i+p-1}}\mathcal{L}_{1}\psi_{j+p-1}\,dx,\qquad
i,j=1,\ldots,n,$ (2.12)

where $i$ counts the rows and $j$ does the columns in the matrix
$\mathrm{M}_{1}$; since the operator $\mathcal{L}_{1}$ acts from
$W_{2}^{2}(\Omega,e^{-\vartheta|x_{d}|}dx)$ into
$L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)$, the functions
$\mathcal{L}_{1}\psi_{j+p-1}$ belong to $L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)$
and this obviously ensures the convergence of the integrals in the above
identity.

By $\mu_{i}$, $i=1,\ldots,N$, we denote different eigenvalues of the matrix
$\mathrm{M}_{1}$ of multiplicities $q_{1},\ldots,q_{N}$. It is clear that
$N\leqslant n$ and $q_{1}+\ldots+q_{N}=n$.

######  Theorem 2.2.

Fix $p\in\\{1,\ldots,m\\}$, $\tau\in\\{-1,+1\\}$ and let
$\Lambda_{p}=\ldots=\Lambda_{p+n-1}$ be an $n$-multiple eigenvalue of the
operator $\mathcal{H}^{\prime}$, where $n\geqslant 1$. There are exactly $N$
poles, counting their orders, of the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ converging to zero as $\varepsilon\to+0$.
These poles, denoted by $k_{ij}(\varepsilon)$, have the asymptotic behavior

$k_{ij}(\varepsilon)=\varepsilon\mu_{i}+O\big{(}\varepsilon^{1+\frac{1}{q_{i}}}\big{)},\qquad
i=1,\ldots,N,\quad j=1,\ldots,q_{i}.$ (2.13)

Asymptotic expansion (2.13) for the poles $k_{ij}$ can be specified in more
details and this will be done in terms of one more matrix
$\mathrm{M}_{2,\tau}$ with entries

$M_{2,\tau}^{ij}:=\frac{1}{2}\int\limits_{\Omega}\overline{\psi_{i+p-1}}(\mathcal{L}_{2}-\mathcal{L}_{1}\mathcal{G}_{p,\tau}\mathcal{L}_{1})\psi_{j+p-1}\,dx\qquad
i,j=1,\ldots,n,$ (2.14)

where $i$ counts the rows and $j$ does the columns in the matrix
$\mathrm{M}_{2,\tau}$. Due to the definition of the operators
$\mathcal{L}_{1}$, $\mathcal{L}_{2}$ , the second term in the integrand in
(2.14) belongs to $L_{2}(\mathds{R},e^{\vartheta|x_{d}|}dx_{d})$ and this
ensures the convergence of the integral. We denote

$Q_{i,\tau}(z):=\frac{\partial\
}{\partial\varepsilon}\det\big{(}z\mathrm{E}-\mathrm{M}_{1}+\varepsilon\mathrm{M}_{2,\tau}\big{)}\bigg{|}_{\varepsilon=0}.$
(2.15)

We stress that if $\Lambda_{p}=\Lambda_{1}$, the matrix $M_{2,\tau}$ and the
function $Q_{i,\tau}$ become independent of $\tau$.

######  Theorem 2.3.

Under the assumptions of Theorem 2.2, we fix $i\in\\{1,\ldots,N\\}$. If
$Q_{i,\tau}(z)$ vanishes identically, then

$k_{ij}(\varepsilon)=\varepsilon\mu_{i}+O\big{(}\varepsilon^{1+\frac{2}{q_{i}}}\big{)},\qquad
i=1,\ldots,N,\quad j=1,\ldots,q_{i}.$ (2.16)

If $Q_{i,\tau}$ is not identically zero, then there exists a fixed non-
negative integer $r_{i,\tau}<q_{i}$ such that

$\gamma_{i,\tau}:=\frac{r_{i,\tau}!}{\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{N}(\mu_{i}-\mu_{j})^{q_{j}}}\frac{d^{r_{i,\tau}}Q_{i,\tau}}{dz^{r_{i,\tau}}}(\mu_{i})\neq
0.$ (2.17)

If $2r_{i,\tau}\geqslant q_{i}$, then

$k_{ij}(\varepsilon)=\varepsilon\mu_{i}+O\big{(}\varepsilon^{1+\frac{1}{r_{i,\tau}}}\big{)},\qquad
i=1,\ldots,N,\quad j=1,\ldots,q_{i}.$ (2.18)

If $2r_{i,\tau}\leqslant q_{i}-1$, then exactly $r_{i,\tau}$ poles $k_{ij}$,
$j=1,\ldots,r_{i,\tau}$ have the asymptotic behavior

$k_{ij}(\varepsilon)=\varepsilon\mu_{i}+O\big{(}\varepsilon^{1+\frac{1}{r_{i,\tau}}}\big{)},\qquad
i=1,\ldots,N,\quad j=1,\ldots,r_{i,\tau},$ (2.19)

while other poles $k_{ij}$, $j=r_{i,\tau}+1,\ldots,q_{i}$, have the asymptotic
behavior

$k_{ij}(\varepsilon)=\varepsilon\mu_{i}+\varepsilon^{1+\frac{1}{q_{i}-r_{i,\tau}}}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}+O\big{(}\varepsilon^{1+\frac{2}{q_{i}-r_{i,\tau}}}\big{)},$
(2.20)

where the branch of the fractional power $z^{\frac{1}{q_{i}-r_{i,\tau}}}$ is
fixed by the condition $1^{\frac{1}{q_{i}-r_{i,\tau}}}=1$ with the branch cut
along the negative real semi-axis.

We give some definitions before we formulate our next result. A pole
$k_{\varepsilon}\in B_{\delta}$ of an operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ corresponds to an eigenvalue
$\Lambda_{p}-k_{\varepsilon}^{2}$ of an operator $\mathcal{H}_{\varepsilon}$
if an associated nontrivial solution $\psi_{\varepsilon}$ to (2.3), (2.7)
belongs to $W_{2}^{2}(\Omega)$. Otherwise it corresponds to a resonance
$\Lambda_{p}-k_{\varepsilon}^{2}$.

Our further results provide conditions allowing to determine whether
resonances or eigenvalues are associated with the poles described in two
previous theorems. We first present the main result on the poles emerging from
the bottom of the essential spectrum.

######  Theorem 2.4.

Let $p=1$ and make the assumptions of Theorems 2.2, 2.3. If
$\operatorname{Re}\mu_{i}>0$, then the poles $k_{ij}$, $j=1,\ldots,q_{i}$
correspond to the eigenvalues
$\lambda_{ij}(\varepsilon)=\Lambda_{p}-k_{ij}^{2}(\varepsilon)$ with the
asymptotic behavior

$\lambda_{ij}(\varepsilon)=\Lambda_{p}-\varepsilon^{2}\mu_{i}^{2}+O\big{(}\varepsilon^{2+\frac{1}{\alpha_{i}}}\big{)}$
(2.21)

with $j=1,\ldots,q_{i}$, where

$\alpha_{i}:=\left\\{\begin{aligned}
&\frac{q_{i}}{2}\qquad\hphantom{\tau}\text{if}\quad Q_{i,\tau}\ \text{vanishes
identically},\\\ &\,r_{i,\tau}\qquad\text{if}\quad 2r_{i,\tau}\geqslant
q_{i}.\end{aligned}\right.$ (2.22)

If $2r_{i,\tau}\leqslant q_{i}-1$, then the eigenvalues $\lambda_{ij}$ still
have asymptotic behavior (2.21) with $\alpha_{i}=r_{i,\tau}$ for
$j=1,\ldots,r_{i,\tau}$, while the asymptotic behaviors for the other
eigenvalues read as

$\lambda_{ij}(\varepsilon)=\Lambda_{p}-\varepsilon^{2}\mu_{i}^{2}-2\varepsilon^{2+\frac{1}{q_{i}-r_{i,\tau}}}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}+O(\varepsilon^{2+\frac{2}{q_{i}-r_{i,\tau}}}),\qquad
j=r_{i,\tau}+1,\ldots,q_{i}.$ (2.23)

If $\operatorname{Re}\mu_{i}<0$, then the poles $k_{ij}$, $j=1,\ldots,q_{i}$
correspond to the resonances
$\lambda_{ij}(\varepsilon)=\Lambda_{p}-k_{ij}^{2}(\varepsilon)$ with
asymptotic expansions (2.21), (2.22), (2.23).

Let $\operatorname{Re}\mu_{i}=0$, $Q_{i,\tau}$ be not identically zero and
$2r_{i,\tau}\leqslant q_{i}-1$. As $j=r_{i,\tau}+1,\ldots,q_{i}$, if

$\operatorname{Re}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}>0,$
(2.24)

then the pole $k_{ij}$ corresponds to an eigenvalue, while if

$\operatorname{Re}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}<0,$
(2.25)

the pole $k_{ij}$ corresponds to a resonance. The asymptotic expansion for
this eigenvalue/resonance is given by (2.23) if $\mu_{i}\neq 0$ and in the
case $\mu_{i}=0$ it reads as

$\lambda_{ij}(\varepsilon)=\Lambda_{p}-\varepsilon^{2+\frac{2}{q_{i}-r_{i,\tau}}}(-\gamma_{i,\tau})^{\frac{2}{q_{i}-r_{i,\tau}}}e^{\frac{4\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}+O\big{(}\varepsilon^{2+\frac{3}{q_{i}-r_{i,\tau}}}\big{)}.$
(2.26)

The next theorem concerns the poles emerging from the internal thresholds in
the essential spectrum. Given $p>1$ such that $\Lambda_{p}>\Lambda_{1}$, and
$\tau\in\\{-1,+1\\}$, by $k_{ij,\tau}=k_{ij,\tau}(\varepsilon)$ we redenote
the corresponding poles $k_{ij}$ of the operator
$\mathcal{R}_{\varepsilon,\tau}$ described in Theorems 2.2, 2.3.

######  Theorem 2.5.

Let $\Lambda_{p}>\Lambda_{1}$ and fix $i\in\\{1,\ldots,N\\}$,
$j\in\\{1,\ldots,q_{i}\\}$, $\tau\in\\{-1,+1\\}$. Let

$\displaystyle\operatorname{Re}\mu_{i}>0$ (2.27) or
$\displaystyle\begin{gathered}\operatorname{Re}\mu_{i}=0,\qquad
Q_{i,\tau}\not\equiv 0,\qquad 2r_{i,\tau}\leqslant q_{i}-1,\\\
j\in\\{r_{i}+1,\ldots,q_{i}\\},\qquad\operatorname{Re}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}>0\end{gathered}$
(2.30)

and

$\displaystyle\tau\operatorname{Im}\mu_{i}<0$ (2.31) or
$\displaystyle\begin{gathered}\operatorname{Im}\mu_{i}=0,\qquad
Q_{i,\tau}\not\equiv 0,\qquad 2r_{i,\tau}\leqslant q_{i}-1,\\\
j\in\\{r_{i}+1,\ldots,q_{i}\\},\qquad\tau\operatorname{Im}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}<0.\end{gathered}$
(2.34)

Then the pole $k_{ij,\tau}(\varepsilon)$ corresponds to an eigenvalue
$\lambda_{ij,\tau}(\varepsilon)=\Lambda_{p}-k_{ij,\tau}^{2}(\varepsilon)$ with
asymptotic expansions (2.21), (2.22), (2.23) if $\mu_{i}\neq 0$ and asymptotic
expansion (2.26) if $\mu_{i}=0$.

Let

$\displaystyle\operatorname{Re}\mu_{i}<0$ (2.35) or
$\displaystyle\begin{gathered}\operatorname{Re}\mu_{i}=0,\qquad
Q_{i,\tau}\not\equiv 0,\qquad 2r_{i,\tau}\leqslant q_{i}-1,\\\
j\in\\{r_{i}+1,\ldots,q_{i}\\},\qquad\operatorname{Re}(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{2\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}<0.\end{gathered}$
(2.38)

Then the pole $k_{ij,\tau}(\varepsilon)$ corresponds to a resonance
$\lambda_{ij,\tau}(\varepsilon)=\Lambda_{p}-k_{ij,\tau}^{2}(\varepsilon)$ with
asymptotic expansions (2.21), (2.22), (2.23) if $\mu_{i}\neq 0$ and asymptotic
expansion (2.26) if $\mu_{i}=0$.

Let $\mu_{i}$ be a simple eigenvalue of the matrix $\mathrm{M}_{1}$ with an
associated eigenvector
$\mathrm{e}_{i}:=(\mathrm{e}_{i,1},\ldots,\mathrm{e}_{i,n})$, $j=1$, condition
(2.27) or (2.30) hold, and

$\displaystyle\tau\operatorname{Im}\mu_{i}<0$ (2.39) or
$\displaystyle\begin{gathered}\operatorname{Im}\mu_{i}=0,\qquad
Q_{i,\tau}\not\equiv 0,\qquad
r_{i,\tau}=0,\qquad\tau\operatorname{Im}\gamma_{i,\tau}<0,\end{gathered}$
(2.41)

and let there exist $s\in\\{1,\ldots,p-1\\}$ such that

$\sum\limits_{t=1}^{n}\int\limits_{\Omega}e^{-K_{t}(0)x_{d}}\overline{\psi_{s}(x^{\prime})}\mathcal{L}_{1}\mathrm{e}_{i,t}\psi_{t-p+1}\,dx\neq
0\quad\text{or}\quad\sum\limits_{t=1}^{n}\int\limits_{\Omega}e^{K_{t}(0)x_{d}}\overline{\psi_{s}(x^{\prime})}\mathcal{L}_{1}\mathrm{e}_{i,t}\psi_{t-p+1}\,dx\neq
0.$ (2.42)

Then the pole $k_{i1,\tau}{(\varepsilon)}$ corresponds to a resonance
$\lambda_{i,\tau}{(\varepsilon)}=\Lambda_{p}-k_{i1,\tau}^{2}(\varepsilon)$
with the asymptotic behavior

$\lambda_{i,\tau}(\varepsilon)=\Lambda_{p}-\varepsilon^{2}\mu_{i}^{2}+O(\varepsilon^{4})$

if $Q_{i,\tau}$ vanishes identically, and

$\lambda_{i,\tau}(\varepsilon)=\Lambda_{p}-\varepsilon^{2}(\mu_{i}-\varepsilon\gamma_{i,\tau})^{2}+O\big{(}|\mu_{i}|\varepsilon^{4}+\varepsilon^{5}\big{)}$

otherwise.

###  2.3 Discussion of the results

In this subsection we discuss the main results formulated in Theorems 2.1,
2.2, 2.3, 2.4, 2.5. The first of them, Theorem 2.1, describes a meromorphic
continuation of the resolvent of the perturbed operator. This continuation is
local and is constructed in the vicinity of the points $\Lambda_{p}$,
$p=1,\ldots,m$. The point $\Lambda_{1}$ is the bottom of the essential
spectrum, see Lemma 2.1 and in vicinity of this point just one meromorphic
continuation is possible. It is introduced as a solution to problem (2.3) with
a specified behaviour at infinity, see (2.4), in terms of an auxiliary
spectral parameter $k$. The right hand side in the equation in (2.3) is not in
the class of compactly supported functions as it is usually assumed for
meromorphic continuations, but an element of a wider space
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|\,dx})$. Here the presence of the weight
$e^{{\vartheta}|x_{d}|}$ means that the elements of latter space in certain
sense decays exponentially as $x_{d}\to\pm\infty$, namely they are represented
as $f=e^{-\frac{{\vartheta}|x_{d}|}{2}}\tilde{f}$, where $\tilde{f}\in
L_{2}(\Omega)$. The final operator providing the meromorphic continuation is
$\mathcal{R}_{\varepsilon,\tau}$ and for $\Lambda_{1}$ it is independent of
$\tau$.

In the vicinity of internal thresholds $\Lambda_{p}>\Lambda_{1}$ in the
essential spectrum, there are _two different_ meromorphic continuations given
by the operators $\mathcal{R}_{\varepsilon,-1}$ and
$\mathcal{R}_{\varepsilon,+1}$. The former describes a meromorphic
continuation from the lower complex half-plane into the upper one, while the
latter does from the upper half-plane into the lower one. In the theory of
self-adjoint operators, usually only the latter continuation from the upper
half-plane into the lower one is studied since it is physically meaning and it
arises while considering a corresponding Cauchy problem for an evolutionary
Schrödinger equation. However, since our perturbing operator is not assumed to
be symmetric, the operator is not necessary self-adjoint. As a result, it can
possess complex eigenvalues in the vicinity of the threshold $\Lambda_{p}$.
These eigenvalues are the poles of the resolvent of the perturbed operator.
And as we shall see below, once we continue meromorphically the resolvent from
the upper half-plane into the lower one, the eigenvalues in the lower half-
plane can become ‘‘invisible’’ for the continuation in the sense that this
continuation has no poles at such eigenvalues. A similar situation can hold
once we continue analytically the resolvent from the lower half-plane into the
upper one. A clear explanation of this phenomenon is due to representations
(2.4), (2.5), (2.6) and (2.7), (2.8), (2.9). Namely, as
$\Lambda_{p}>\Lambda_{1}$, the functions $\phi_{\varepsilon,j}^{\pm}$,
$j=1,\ldots,p-1$, in (2.8) behave at infinity as
$\phi_{\varepsilon,j}^{\pm}(x_{d})\sim e^{-\tau K_{j}(k)|x_{d}|}$. In view of
obvious identities

$K_{j}(k)=-\mathrm{i}\sqrt{\Lambda_{p}-\Lambda_{j}}+\frac{\mathrm{i}}{2\sqrt{\Lambda_{p}-\Lambda_{j}}}k^{2}+O(k^{4}),\qquad
k\to 0,$ (2.43)

the exponents $e^{-\tau K_{j}(k)|x_{d}|}$ decay only if
$\tau\operatorname{Im}k^{2}<0$. Depending on $\tau$, the latter condition
means that in general only the eigenvalues either in the upper or lower
complex half-plane can serve as poles of the meromorphic continuation of the
resolvent of the operator $\mathcal{H}_{\varepsilon}$. This is a main reason
why we deal with both meromorphic continuations, in contrast to the case of
self-adjoint operators with symmetric perturbations.

Theorems 2.2, 2.3 describe the poles of the meromorphic continuations of the
resolvent in the vicinity of the thresholds $\Lambda_{p}$ in the essential
spectrum. The first theorem states that in the vicinity of an $n$-multiple
threshold $\Lambda_{p}=\ldots=\Lambda_{p+n-1}$ there exist exactly $n$ poles
of the operator $\mathcal{R}_{\varepsilon,\tau}$ counting their orders. We
stress that here we count the _orders_ of the poles and not their
_multiplicities_ , that is, not _the number of associated linear independent
solutions_ to problem (2.3) with $f=0$, $k=k_{ij}(\varepsilon)$. The
multiplicity of each pole does not exceed its order; this can be shown by the
technique used in the proofs of Lemmata 6.2, 6.3 in [6] and Lemmata 6.2, 6.3
in [7]. However, in general, the multiplicities and the orders coincide only
if the perturbation $\mathcal{L}(\varepsilon)$ is symmetric. The reason is
that in the general case of a non-symmetric perturbation, in a certain matrix
controlling the structure of the poles $k_{ij}$, a non-diagonal Jordan block
can arise and this gives rise to adjoint vectors instead of the eigenvectors,
see Section 5 and the calculations involving matrix
$\mathrm{M}_{\varepsilon,\tau}$. Of course, the multiplicity of each pole
$k_{ij,\tau}$ is at least one. In particular, if all poles $k_{ij,\tau}$ are
different for a fixed $\tau$, the total multiplicity is equal to $n$. Theorem
2.2 provides leading terms in the asymptotic expansions for the poles
$k_{ij,\tau}$, while Theorem 2.3 specifies these expansions. In some cases it
just improves the estimate for the error terms, see (2.16), (2.18), (2.19),
while in some cases, a next-to-leading term in the expansions can be found,
see (2.20). Theorems 2.2, 2.3 treat a general case, when the eigenvalues of
the matrix $\mathrm{M}_{1}$ are of arbitrary multiplicities and no extra
assumptions are made for the matrix $\mathrm{M}_{2,\tau}$. In an important
particular case, when $\mu_{i}$ is a simple eigenvalue of the matrix
$\mathrm{M}_{1}$, we have $q_{i}=1$ and $r_{i,\tau}=0$. In this case there
exists just one pole $k_{i1,\tau}$ with asymptotic behavior (2.13) and
expansion (2.16), (2.20) can be applied, which yields that

$k_{i1,\tau}(\varepsilon)=\varepsilon\mu_{i}-\varepsilon^{2}\gamma_{i,\tau}+O(\varepsilon^{3}).$

If $n=1$, that is, $\Lambda_{p}$ is a simple eigenvalue of the operator
$\mathcal{H}^{\prime}$, the above expansions can be specified as follows:

$k_{i1,\tau}(\varepsilon)=-\frac{\varepsilon}{2}\int\limits_{\Omega}\overline{\psi_{p}}\mathcal{L}_{1}\psi_{p}\,dx-\frac{\varepsilon^{2}}{2}\int\limits_{\Omega}\overline{\psi_{p}}(\mathcal{L}_{2}-\mathcal{L}_{1}\mathcal{G}_{p,\tau}\mathcal{L}_{1})\psi_{p}\,dx+O(\varepsilon^{3}).$
(2.44)

We also observe that since the operator $\mathcal{R}_{\varepsilon,\tau}$ is
independent of $\tau$ if $p=1$, in the general situation there are only $n$
poles in the vicinity of the bottom $\Lambda_{1}$ of the essential spectrum.
In the vicinity of internal thresholds $\Lambda_{p}>\Lambda_{1}$, the
operators $\mathcal{R}_{\varepsilon,\tau}$ depend on $\tau$ and this is why
there are $2n$ poles in the vicinity of $\Lambda_{p}$. In particular, if
$\Lambda_{p}$ is an $n$-multiple eigenvalue of the operator
$\mathcal{H}^{\prime}$, there can be $2n$ different simple eigenvalues of the
operator $\mathcal{H}_{\varepsilon}$ converging to $\Lambda_{p}$, see examples
in Subsection 3.3.

The above discussed poles of the operators $\mathcal{R}_{\varepsilon,\tau}$
correspond either to the eigenvalues or resonances depending on the behavior
of the associated non-trivial solutions. This behaviour is completely
described by formulae (2.7), (2.8) and we just need to identify whether the
function $\psi_{\varepsilon}$ decays exponentially at infinity or not. In the
former case we deal with an eigenvalue, otherwise with a resonance. As we see,
the functions $\phi_{\varepsilon,j}^{\pm}(x_{d})$ decay exponentially as
$j\geqslant p+n$ no matter how the corresponding pole looks like. But for
$j=p,\ldots,p+n-1$ these functions behave at infinity as
$\phi_{\varepsilon,j}^{\pm}(x_{d},k)\sim e^{-k_{\varepsilon}|x_{d}|}$. These
functions decay exponentially as $\operatorname{Re}k_{\varepsilon}>0$, is
periodic as $\operatorname{Re}k_{\varepsilon}=0$ and grows exponentially as
$\operatorname{Re}k_{\varepsilon}<0$. If $\Lambda_{p}>\Lambda_{1}$, we also
have to control the behavior of the functions
$\phi_{\varepsilon,j}^{\pm}(x_{d})$ with $j=1,\ldots,p-1$. This is easily done
by identities (2.43): the functions $\phi_{\varepsilon,j}^{\pm}(x_{d})$,
$j=1,\ldots,p-1$, decay exponentially if
$\tau\operatorname{Im}k_{\varepsilon}^{2}<0$, are periodic if
$\operatorname{Im}k_{\varepsilon}^{2}=0$ and grow exponentially if
$\tau\operatorname{Im}k_{\varepsilon}^{2}>0$. All discussed conditions can be
checked by means of asymptotic expansions provided by Theorems 2.2, 2.3 for a
given pole. And exactly this is done in the proof of Theorems 2.4, 2.5.
Conditions in Theorem 2.4 are aimed at checking the sign of the real part of a
given pole and proving at the same time that at least one of the coefficients
$c_{\varepsilon,j}^{\pm}$, $j=1,\ldots,n$, in (2.8) is non-zero. Similar
conditions (2.27), (2.30), (2.31), (2.34), (2.35), (2.38) ensure that the
functions $\phi_{\varepsilon,j}^{\pm}$, $j=1,\ldots,p-1$ decay exponentially,
that is, $\tau\operatorname{Im}k_{\varepsilon}^{2}<0$, while for
$j=p,\ldots,p+n-1$, these functions demonstrate either an exponential decay or
an exponential growth. Conditions (2.39), (2.41), (2.42) describe a more
gentle situation. Namely, here the real part of the pole is negative and the
functions $\phi_{\varepsilon,j}^{\pm}$, $j=p,\ldots,p+n-1$, decay
exponentially. However, $\tau\operatorname{Im}k_{\varepsilon}^{2}>0$ and this
means that the functions $\phi_{\varepsilon,j}^{\pm}$, $j=1,\ldots,p-1$, can
grow exponentially. This is true, once we guarantee that at least one of the
coefficients $c_{\varepsilon,j}^{\pm}$, $j=1,\ldots,p-1$, is non-zero. This is
indeed the case thanks to condition (2.42). The asymptotic expansions for the
eigenvalues and the resonances provided in Theorems 2.4, 2.5 are implied
immediately by the formula
$\lambda_{\varepsilon}=\Lambda_{p}-k_{\varepsilon}^{2}$ relating the
eigenvalues/resonances with a pole $k_{\varepsilon}$ and the asymptotic
expansions for the poles stated in Theorems 2.2, 2.3.

Let us briefly discuss the main ideas underlying our main results. First, we
rather straightforwardly construct the meromorphic continuation for the
resolvent of the unperturbed operator. Namely, we find explicitly the
projection of the solution to problem (2.3) on the eigenfunctions $\psi_{j}$,
$j=1,\ldots,m$, and study then the properties of the coefficients in this
projection and of the remaining orthogonal part in the solution. Once such
continuation is constructed, for proving our main results, we apply an
approach being a modification of the technique suggested in [5], [8], [25],
[26]. The idea is to regard the perturbation as a right-hand side and to apply
then the meromorphic continuation of the unperturbed operator. After some
simple calculations this leads us to an operator equation with a certain
finite rank perturbation. Resolving this equation, we rather easily succeed to
construct the meromorphic continuation for the perturbed operator and identify
its poles as solutions to a nonlinear eigenvalue problem for some explicitly
calculated matrix depending also on the small parameter. Analysing then this
problem by means of methods from the theory of complex functions, we study the
existence of the poles and their asymptotic behavior.

Although we restrict ourselves by considering second order differential
operators, our approach can be also adapted for certain operators of higher
order. However, general higher order operators can have a richer spectral
structure of the edges in the essential spectrum and there can be more
complicated scenarios of their bifurcations under perturbations, see, for
instance, [55]. This is why the case of second order operators deserves a
separate study, what is done in the present work. Our results are of general
nature and are applicable to wide classes of unperturbed operators and
perturbations. In the next section we discuss some possible examples of both
unperturbed operators and perturbations as well as some specific examples
motivated by physical models.

##  3 Examples

In this section we provide examples demonstrating our main results.

###  3.1 Unperturbed operator

Here we discuss some examples of the unperturbed operator, namely, of the
operator $\mathcal{H}$. This is a general self-adjoint second order
differential operator and it includes such classical operators as a
Schrödinger operator:

$\mathcal{H}=-\Delta+A_{0},\qquad A_{0}=A_{0}(x^{\prime}),$

a magnetic Schrödinger operator:

$\mathcal{H}=(\mathrm{i}\nabla_{x^{\prime}}+{A})^{2}-\frac{\partial^{2}\
}{\partial
x_{d}^{2}}+A_{0},\qquad{A=(A_{1},\ldots,A_{d-1})},\qquad{A_{j}=A_{j}(x^{\prime})},\qquad
A_{0}=A_{0}(x^{\prime}),$

a Schrödinger operator with metric:

$\mathcal{H}=-\sum\limits_{i,j=1}^{d-1}\frac{\partial\ }{\partial
x_{i}}A_{ij}\frac{\partial\ }{\partial x_{j}}-\frac{\partial^{2}\ }{\partial
x_{d}^{2}}+A_{0},\qquad A_{ij}=A_{ij}(x^{\prime}),\qquad
A_{0}=A_{0}(x^{\prime}).$

All these operators are considered in a tubular domain
$\Omega=\omega\times\mathds{R}$. If $\omega=\mathds{R}^{d-1}$, then the domain
$\Omega$ becomes an entire space $\mathds{R}^{d}$. If $\omega$ is a bounded
domain, not necessary connected, then $\Omega$ is an infinite cylinder, which
is to be regarded as a quantum waveguide if the Dirichlet condition is imposed
on its boundary and as an acoustic waveguide if the boundary is subject to the
Neumann condition. Further examples of unbounded domains $\omega$ are also
possible. For instance, if $\omega$ is the half-space
$\omega:=\\{x^{\prime}:\,x_{j}>0\\}$ in $\mathds{R}^{d-1}$ for some
$j=1,\ldots,d-1$, the domain $\Omega$ becomes the half-space
$\\{x:\,x_{j}>0\\}$ in $\mathds{R}^{d}$. We can also consider a more
complicated domain $\omega:=\\{x^{\prime}:x_{d-1}<h(x_{1},\ldots,x_{d-2})\\}$
for some smooth function $h$, then
$\Omega=\\{x:\,x_{d-1}<h(x_{1},\ldots,x_{d-2}),x_{d}\in\mathds{R}\\}$.

###  3.2 Perturbation

In this subsection we discuss possible examples of the perturbing operator
$\mathcal{L}(\varepsilon)$. The first example is a second order differential
operator:

$\mathcal{L}(\varepsilon)=\sum\limits_{i,j=1}^{n}{\Upsilon}_{ij}(x,\varepsilon)\frac{\partial^{2}\
}{\partial x_{i}\partial
x_{j}}+\sum\limits_{j=1}^{n}{\Upsilon}_{j}(x,\varepsilon)\frac{\partial\
}{\partial x_{j}}+{\Upsilon}_{0}(x,\varepsilon).$ (3.1)

Here ${\Upsilon}_{ij},\,{\Upsilon}_{j},\,{\Upsilon}_{0}\in L_{\infty}(\Omega)$
are some functions, not necessarily real-valued, satisfying the
representations

${\Upsilon}_{\natural}(x,\varepsilon)={\Upsilon}_{\natural}^{(1)}(x)+\varepsilon{\Upsilon}_{\natural}^{(2)}(x)+\varepsilon^{2}{\Upsilon}_{\natural}^{(3)}(x,\varepsilon),\qquad\natural=ij,j,0,$
(3.2)

where ${\Upsilon}_{\natural}^{(s)}\in L_{\infty}(\Omega)$ are some functions
obeying the estimates:

$\|{\Upsilon}_{\natural}^{(s)}e^{2{\vartheta}|x_{d}|}\|_{L_{\infty}(\Omega)}<C,$

and $\vartheta$ and $C$ are some fixed positive constant independent of
$\varepsilon$. In this case the operators $\mathcal{L}_{s}$ read as

$\mathcal{L}_{s}:=\sum\limits_{i,j=1}^{n}{\Upsilon}_{ij}^{(s)}\frac{\partial^{2}\
}{\partial x_{i}\partial
x_{j}}+\sum\limits_{j=1}^{n}{\Upsilon}_{j}^{(s)}\frac{\partial\ }{\partial
x_{j}}+{\Upsilon}_{0}^{(s)},\qquad s=1,2,3.$

Particular cases of this example are small potential, small magnetic field,
small metric.

The second example is an integral operator of the form

$(\mathcal{L}(\varepsilon)u)(x,\varepsilon)=\int\limits_{\Omega}{J}(x,y,\varepsilon)u(y)\,dy,$

where ${J}\in L_{2}(\Omega\times\Omega)$ is some kernel, not necessarily real-
valued and symmetric and satisfying the representation:

${J}(x,y,\varepsilon)={J}_{1}(x,y)+\varepsilon{J}_{2}(x,y)+\varepsilon^{2}{J}_{3}(x,y,\varepsilon),$

where $L_{i}$ are some functions obeying the estimates:

$\int\limits_{\Omega\times\Omega}|{J}_{i}|^{2}e^{{vt}(|x_{d}|+|y_{d}|)}\,dx\,dy<C,\qquad
i=1,2,3,$

and ${vt}$ and $C$ are some fixed positive constant independent of
$\varepsilon$.

The third example is a localized $\delta$-interaction with a complex-valued
density. Namely, let $S\subset\Omega$ be a manifold of codimension $1$ and of
smoothness $C^{3}$. We assume that it is compact and has no edge. The
perturbed operator in question is

$\mathcal{H}_{\varepsilon}=-\Delta+\varepsilon\beta\delta(x-S),$

which acts as $\mathcal{H}_{\varepsilon}u=-\Delta u$ on the domain
$\mathfrak{D}(\mathcal{H}_{\varepsilon})$ formed by the functions $u\in
W_{2}^{2}(\Omega\setminus S)\cap W_{2}^{1}(\Omega)$ obeying the boundary
condition $\mathcal{B}u=0$ on $\partial\Omega$ and the boundary conditions

$[u]_{S}=0,\qquad\left[\frac{\partial
u}{\partial\nu}\right]_{S}=\varepsilon\beta u.$

Here $[u]_{S}$ denotes the jump of the function on $S$, namely,

$[u]_{S}:={\lim\limits_{t\to 0+}\big{(}u(\cdot+t\nu)-u(\cdot-t\nu)\big{)},}$

and $\nu$ is the unit normal to $S$ directed outside the domain enveloped by
$S$. By $\beta$ we denote some complex-valued function defined on $S$,
uniformly bounded and belonging to $C^{2}(S)$ . Such operator does not satisfy
our assumptions for $\mathcal{L}(\varepsilon)$ since now the perturbation
changes the domain. However, it is possible to reduce the perturbed operator
to another one obeying needed assumptions and having the same eigenvalues and
resonances. Namely, thanks to the made assumptions on the manifold $S$, in a
small vicinity of the $S$ we can introduce a new variable $\rho$ being the
distance from a point to $S$ measured along the normal $\nu$. This variable is
well-defined at least in the neighbourhood
$\\{x:\,\operatorname{dist}(x,S)<\rho_{0}\\}$ of $S$, where $\rho_{0}$ is some
fixed number. Let $\chi=\chi(\rho)$ be an infinitely differentiable function
such that $\chi(\rho)=\frac{|\rho|}{2}$ as $|\rho|<\frac{\rho_{0}}{3}$ and
$\chi(\rho)=0$ as $|\rho|>\frac{2\rho_{0}}{3}$. By $\mathcal{U}_{\varepsilon}$
we denote the multiplication operator
$\mathcal{U}_{\varepsilon}u:=(1+\varepsilon\beta\chi)^{-1}u$. It is
straightforward to check that this operator maps the domain of the operator
$\mathcal{H}_{\varepsilon}$ onto the space $\\{u\in
W_{2}^{2}(\Omega):\,\mathcal{B}u=0\ \text{on}\ \partial\Omega\\}$. This space
serves as the domain for an operator
$\tilde{\mathcal{H}}_{\varepsilon}:=\mathcal{U}_{\varepsilon}\mathcal{H}_{\varepsilon}\mathcal{U}_{\varepsilon}^{-1}$.
It is straightforward to confirm that the differential operator for the latter
operator reads as

$\tilde{\mathcal{H}}_{\varepsilon}=-\Delta+\varepsilon\mathcal{L}(\varepsilon),\qquad\mathcal{L}(\varepsilon):=-2(1+\varepsilon\beta\chi)^{-1}\nabla\beta\chi\cdot\nabla-(1+\varepsilon\beta\chi)^{-1}\Delta\beta\chi$

and we see that a first order differential operator $\mathcal{L}(\varepsilon)$
is a particular case of operator (3.1). It is also clear that the operators
$\mathcal{H}_{\varepsilon}$ and $\tilde{\mathcal{H}}_{\varepsilon}$ have the
same eigenvalues and resonances since the operator $\mathcal{U}_{\varepsilon}$
does not change the behavior of the functions at infinity. Hence, we can study
the eigenvalues and resonances of the operator
$\tilde{\mathcal{H}}_{\varepsilon}$ and transfer then the results to the
operator $\mathcal{H}_{\varepsilon}$.

Our fourth example is a geometric perturbation. Namely, let $\omega$ have a
non-empty boundary, then the same is true for $\Omega$. By $\Gamma$ we denote
a bounded subset of the boundary $\partial\Omega$. Let $\rho$ be a distance to
a point measured along the outward normal to $\partial\Omega$ and $h\in
C^{2}(\partial\Omega)$ be some real function defined on $\Gamma$ and compactly
supported in $\Gamma$. Then we consider a domain $\Omega_{\varepsilon}$
obtained by a small variation of the part $\Gamma$ of the boundary
$\partial\Omega$. Namely, $\Omega_{\varepsilon}$ is a domain with the
following boundary

$\partial\Omega_{\varepsilon}:=(\partial\Omega\setminus\Gamma)\cup\\{x:\,\rho=\varepsilon
h\\}.$

In such domain we consider an operator with differential expression (2.1)
subject to the Dirichlet boundary condition or Neumann condition. We assume
that all the coefficients in the differential expression depend on
$x^{\prime}$ only and are infinitely differentiable. Such perturbed operator
does not fit our scheme since here the domain $\Omega_{\varepsilon}$ depends
on $\varepsilon$. However, as in the previous example, it is possible to
transform such operator to another one fitting our assumptions. Namely, let
$\chi=\chi(x)$ be an infinitely differentiable cut-off function equalling to
one in some fixed sufficiently small $d$-dimensional neighbourhood of $\Gamma$
and vanishing outside some bigger neighbourhood. In this bigger neighbourhood
we introduce local coordinates $(P,\rho)$, where $P\in\partial\Omega$. A point
$x$ is recovered from $(P,\rho)$ by measuring the distance $\rho$ along the
outward normal to $\partial\Omega$ at the point $P$. Then we define a mapping
$\mathcal{P}$ by the following rule: for each point $x$, we find corresponding
$(P,\rho)$ and the action of the mapping is a point corresponding to
$(P,\rho-\varepsilon h(P))$. We introduce new coordinates by the formula
$\tilde{x}:=x(1-\chi(x))+\varepsilon\chi(x)\mathcal{P}(x)$. It is easy to see
that these coordinates are well-defined provided $\varepsilon$ is small enough
and after passing to these new coordinates, the domain $\Omega_{\varepsilon}$
transforms into $\Omega$, while the operator $\mathcal{H}_{\varepsilon}$
becomes $\mathcal{H}+\mathcal{L}(\varepsilon)$, where
$\mathcal{L}(\varepsilon)$ is some second order differential operator of form
(3.1) with compactly supported coefficients obeying (3.2).

###  3.3 Emerging poles for particular models

In this section we apply Theorems 2.4, 2.5 to some simple two- and three-
dimensional operators motivated by an interesting physical background.

####  3.3.1 Planar waveguide

The first model is an infinite planar waveguide modeled by the Dirichlet
Laplacian. Namely, we let $d=2$, $\omega:=(0,\pi)$, and
$\mathcal{H}^{\prime}=-\frac{d^{2}\ }{dx_{2}^{2}}$ subject to the Dirichlet
boundary condition. Then $\Omega:=\\{x:\,0<x_{1}<\pi\\}$ is an infinite strip
and $\mathcal{H}=-\Delta$ is the Dirichlet Laplacian in $\Omega$. As a
perturbation, we choose a complex-valued potential of the form
$\mathcal{L}(\varepsilon):=V_{1}+\varepsilon V_{2}$, where $V_{i}=V_{i}(x)$
are some continuous compactly supported complex-valued functions. The operator
$\mathcal{H}^{\prime}$ has a purely discrete spectrum formed by simple
eigenvalues $\Lambda_{p}:=p^{2}$, $p\in\mathds{N}$, and the associated
eigenfunctions normalized in $L_{2}(0,\pi)$ are
$\psi_{p}(x_{d}):=\frac{\sqrt{2}}{\sqrt{\pi}}\sin px_{1}$.

This situation models a slab optical waveguide of a finite width, where the
cladding in direction $x_{1}$ imposes zero boundary conditions at $x_{1}=0$
and $x_{1}=\pi$. Assuming that the waveguide is infinite in the second
direction, a paraxial diffraction of an incident beam can be described using
the normalized equation in the form (see e.g. [36])

$\mathrm{i}\partial_{z}\Phi+\Delta\Phi-\varepsilon V_{\varepsilon}(x)\Phi=0,$
(3.3)

where $\Phi(x_{1},x_{2},z)$ corresponds to complex amplitude of the electrical
field, the optical potential $V_{\varepsilon}(x)$ describes a weak localized
modulation of the complex-valued refractive index, and $z$ is the direction of
propagation of the pulse. For stationary modes $\Phi=e^{-\mathrm{i}\lambda
z}\psi$, where $-\lambda$ has the meaning of propagation constant, equation
(3.3) reduces to the eigenvalue problem in the above described planar
waveguide for the equation

$-\Delta\psi+\varepsilon V_{\varepsilon}\psi=\lambda\psi.$

This is exactly the mathematical model we formulated above once we let
$V_{\varepsilon}=V_{1}+\varepsilon V_{2}$. Let us consider the bifurcation of
the thresholds $p^{2}$ under the presence of a small localized potential
$V_{\varepsilon}$.

For $p=1$, there is just one pole $k_{\varepsilon}$ and according formula
(2.44), its asymptotic behavior reads as

$k_{\varepsilon}=-\frac{\varepsilon}{\pi}\int\limits_{\Omega}V_{1}(x)\sin^{2}x_{1}\,dx-\frac{\varepsilon^{2}}{\pi}\int\limits_{\Omega}\big{(}V_{2}(x)\sin^{2}x_{1}-V_{1}(x)U(x)\sin
x_{1}\big{)}\,dx+O(\varepsilon^{3}),$ (3.4)

where $U:=\mathcal{G}_{1,\tau}(V_{1}\sin x_{1})$. This function is given by
formula (2.10). The term
$U^{\bot}:=((\mathcal{H}^{\bot}-\Lambda_{1})^{-1}f^{\bot})(x)$ with
$f=V_{1}\sin x_{1}$ and $f^{\bot}$ defined by (2.11) solves the boundary value
problem

$(-\Delta-1)U^{\bot}=V_{1}\sin x_{1}-\frac{2}{\pi}\sin
x_{1}\int\limits_{0}^{\pi}V_{1}(t_{1},x_{2})\sin^{2}t_{1}\,dt_{1}\quad\text{in}\quad\Omega,\qquad
U^{\bot}=0\quad\text{on}\quad\partial\Omega.$

This problem can be solved explicitly by the separation of variables and this
gives the final formula for $U$:

$\displaystyle U(x)=$ $\displaystyle-\frac{1}{\pi}\sin
x_{1}\int\limits_{\Omega}|x_{2}-t_{2}|V_{1}(t_{1},t_{2})\sin^{2}t_{1}\,dt$
$\displaystyle+\sum\limits_{j=2}^{\infty}\frac{\sin
jx_{1}}{\pi\sqrt{j^{2}-1}}\int\limits_{\Omega}e^{-\sqrt{j^{2}-1}|x_{2}-t_{2}|}V_{1}(t_{1},t_{2})\sin
t_{1}\sin jt_{1}\,dt.$

Hence,

$\displaystyle\int\limits_{\Omega}\big{(}V_{2}(x)$
$\displaystyle\sin^{2}x_{1}-V_{1}(x)U(x)\sin x_{1}\big{)}\,dx$
$\displaystyle=$
$\displaystyle\int\limits_{\Omega}V_{2}(x)\sin^{2}x_{1}\,dx+\frac{1}{\pi}\int\limits_{\Omega^{2}}|x_{2}-t_{2}|V_{1}(x)V_{1}(t)\sin^{2}t_{1}\sin^{2}x_{1}\,dt\,dx$
$\displaystyle-\sum\limits_{j=2}^{\infty}\frac{1}{\pi\sqrt{j^{2}-1}}\int\limits_{\Omega^{2}}e^{-\sqrt{j^{2}-1}|x_{2}-t_{2}|}V_{1}(t)V_{1}(x)\sin
t_{1}\sin jt_{1}\sin x_{1}\sin jx_{1}\,dt\,dx.$

Now we apply Theorem 2.4 and we see that if

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}x_{1}\,dx<0$

or

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}x_{1}\,dx=0\quad\text{and}\quad\operatorname{Re}\int\limits_{\Omega}\big{(}V_{2}(x)\sin^{2}x_{1}-V_{1}(x)U(x)\sin
x_{1}\big{)}\,dx<0,$

then the pole $k_{\varepsilon}$ corresponds to an eigenvalue
$\lambda_{\varepsilon}=1-k_{\varepsilon}^{2}$. And if

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}x_{1}\,dx>0$

or

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}x_{1}\,dx=0\quad\text{and}\quad\operatorname{Re}\int\limits_{\Omega}\big{(}V_{2}(x)\sin^{2}x_{1}-V_{1}(x)U(x)\sin
x_{1}\big{)}\,dx>0,$

then the pole $k_{\varepsilon}$ corresponds to a resonance
$\lambda_{\varepsilon}=1-k_{\varepsilon}^{2}$. The asymptotic expansion for
this eigenvalue/resonance is given by (2.23), (2.26) but it is more
straightforward to find it by (3.4) and the above formula for
$\lambda_{\varepsilon}$.

We proceed to the case $p>1$. Here we again apply formula (2.44) to obtain

$k_{\varepsilon,\tau}=-\frac{\varepsilon}{\pi}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx-\frac{\varepsilon^{2}}{\pi}\int\limits_{\Omega}\big{(}V_{2}(x)\sin^{2}px_{1}-V_{1}(x)U_{\tau}(x)\sin
px_{1}\big{)}\,dx+O(\varepsilon^{3}),$ (3.5)

where $U_{\tau}:=\mathcal{G}_{\tau}V_{1}\sin px_{1}$ is given by formula
(2.10). The term
$U_{\tau}^{\bot}:=((\mathcal{H}^{\bot}-\Lambda_{p})^{-1}f^{\bot})(x)$ with
$f=V_{1}\sin px_{1}$ and $f^{\bot}$ defined by (2.11) solves the boundary
value problem

$(-\Delta-p^{2})U_{\tau}^{\bot}=V_{1}\sin px_{1}-\frac{2}{\pi}\sin
px_{1}\int\limits_{0}^{\pi}V_{1}(t_{1},x_{2})\sin^{2}pt_{1}\,dt_{1}\quad\text{in}\quad\Omega,\qquad
U_{\tau}^{\bot}=0\quad\text{on}\quad\partial\Omega.$

The solution is again given by the separation of variables and a final formula
for $U_{\tau}$ reads as

$\displaystyle U_{\tau}(x)=\sum\limits_{j=1}^{p-1}\frac{\mathrm{i}\tau\sin
jx_{1}}{\pi\sqrt{p^{2}-j^{2}}}\int\limits_{\mathds{R}}e^{\mathrm{i}\tau\sqrt{p^{2}-j^{2}}|x_{2}-t_{2}|}U_{j}(t_{2})\,dt_{2}$
(3.6)
$\displaystyle\hphantom{U_{\tau}(x)=\sum\limits_{j=1}^{p-1}}-\frac{1}{\pi}\sin
px_{1}\int\limits_{\mathds{R}}|x_{2}-t_{2}|U_{p}(t_{2})\,dt_{2}$
$\displaystyle\hphantom{U_{\tau}(x)=\sum\limits_{j=1}^{p-1}}+\sum\limits_{j=p+1}^{\infty}\frac{\sin
jx_{1}}{\pi\sqrt{j^{2}-p^{2}}}\int\limits_{\mathds{R}}e^{-\sqrt{j^{2}-p^{2}}|x_{2}-t_{2}|}U_{j}(t_{2})\,dt_{2}.$
$\displaystyle U_{j}(x_{2}):=\int\limits_{0}^{\pi}V_{1}(t_{1},x_{2})\sin
pt_{1}\sin jt_{1}\,dt,\qquad j\neq p,$

Then we get:

$\displaystyle\int\limits_{\Omega}V_{1}(x)U_{\tau}(x)\sin px_{1}\,dx=$
$\displaystyle\sum\limits_{j=1}^{p-1}\int\limits_{\mathds{R}^{2}}\frac{\mathrm{i}\tau
e^{\mathrm{i}\tau\sqrt{p^{2}-j^{2}}|x_{2}-t_{2}|}}{\pi\sqrt{p^{2}-j^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dx_{2}\,dt_{2}$
(3.7)
$\displaystyle\hphantom{\sum\limits_{j=1}^{p-1}}-\frac{1}{\pi}\int\limits_{\mathds{R}^{2}}|x_{2}-t_{2}|U_{p}(x_{2})U_{p}(t_{2})\,dx_{2}\,dt_{2}$
$\displaystyle+\sum\limits_{j=p+1}^{\infty}\int\limits_{\mathds{R}^{2}}\frac{e^{-\sqrt{j^{2}-p^{2}}|x_{2}-t_{2}|}}{\pi\sqrt{j^{2}-p^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}.$

Now we can apply Theorem 2.5 for $\tau=+1$ and $\tau=-1$ and to determine
whether the poles $k_{\varepsilon,\tau}$ correspond to eigenvalues or
resonances. As we see, in a general situation we can have two eigenvalues or
two resonances or one eigenvalue and one resonance. Let us show that each of
these situations is possible.

First of all we observe that in notations of Theorem 2.5 we have $N=1$,
$q_{i}=1$, $r_{i}=0$,

$\displaystyle\mu_{1}=-\frac{1}{\pi}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx,$
$\displaystyle(-\gamma_{i,\tau})^{\frac{1}{q_{i}-r_{i,\tau}}}e^{\frac{\pi\mathrm{i}}{q_{i}-r_{i,\tau}}(j-r_{i,\tau})}=-\frac{1}{\pi}\int\limits_{\Omega}\big{(}V_{2}(x)\sin^{2}px_{1}-V_{1}(x)U_{1,\tau}(x)\sin
px_{1}\big{)}\,dx.$

Assume now that $V_{1}$ is a complex-valued potential such that

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx>0.$

Then condition (2.35) is satisfied and both poles $k_{\varepsilon,\tau}$,
$\tau=\\{-1,+1\\}$, correspond to resonances.

If

$\operatorname{Re}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx<0,\qquad\operatorname{Im}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx\neq
0,$

then conditions (2.27), (2.31) are satisfied with

$\tau=\operatorname{sgn}\operatorname{Im}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx$

and for such $\tau$, the pole $k_{\varepsilon,\tau}$ corresponds to an
eigenvalue. If in addition,

$\int\limits_{\Omega}e^{-\mathrm{i}\sqrt{p^{2}-s^{2}}x_{2}}V_{1}(x)\sin
sx_{1}\sin px_{1}\,dx\neq
0\quad\text{or}\quad\int\limits_{\Omega}e^{\mathrm{i}\sqrt{p^{2}-s^{2}}x_{2}}V_{1}(x)\sin
sx_{1}\sin px_{1}\,dx\neq 0,$

for some $s\in\\{1,\ldots,p-1\\}$, then conditions (2.39), (2.42) are
satisfied and the pole $k_{\varepsilon,\tau}$ with

$\tau=-\operatorname{sgn}\operatorname{Im}\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx$

corresponds to a resonance.

In order to realize a situation with two eigenvalues, we assume that $V_{2}=0$
and consider a special class of $\mathcal{PT}$-symmetric potentials $V_{1}$.
Namely, we suppose that

$V_{1}(x)=W_{1}(x)+\mathrm{i}W_{2}(x),$ (3.8)

where $W_{1}$, $W_{2}$ are real-valued compactly supported potentials with
certain parity:

$W_{1}(x_{1},-x_{2})=W_{1}(x_{1},x_{2}),\qquad
W_{2}(x_{1},-x_{2})=-W_{2}(x_{1},x_{2}).$ (3.9)

These assumptions yield that the operator $\mathcal{H}_{\varepsilon}$ is
$\mathcal{PT}$-symmetric (or partially $\mathcal{PT}$-symmetric using the
terminology from [68]). They also imply immediately that

$\int\limits_{\Omega}V_{1}(x)\sin^{2}px_{1}\,dx=\int\limits_{\Omega}W_{1}(x)\sin^{2}px_{1}\,dx$

and we assume that

$\int\limits_{\Omega}W_{1}(x)\sin^{2}px_{1}\,dx<0.$ (3.10)

It follows from assumptions (3.9) and the definition of the functions $U_{j}$
in (3.6) that these functions are given by the formulae

$U_{j}(x_{2}):=W_{1,j}(x_{2})+\mathrm{i}W_{2,j}(x_{2}),\qquad
W_{s,j}(x_{2}):=\int\limits_{0}^{\pi}W_{s}(t_{1},x_{2})\sin pt_{1}\sin
jt_{1}\,dt,\qquad s=1,2,$

and the functions $W_{1,j}$ are even, while $W_{2,j}$ are odd. As $j\geqslant
p+1$, by making the change of the variables $x_{2}\mapsto-x_{2}$,
$t_{2}\mapsto-t_{2}$ in the integrals in the second sum in (3.7), we get:

$\int\limits_{\mathds{R}^{2}}\frac{e^{-\sqrt{j^{2}-p^{2}}|x_{2}-t_{2}|}}{\pi\sqrt{j^{2}-p^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}=\int\limits_{\mathds{R}^{2}}\frac{e^{-\sqrt{j^{2}-p^{2}}|x_{2}-t_{2}|}}{\pi\sqrt{j^{2}-p^{2}}}\overline{U_{j}(x_{2})}\overline{U_{j}(t_{2})}\,dt_{2}\,dx_{2}$

and hence,

$\operatorname{Im}\sum\limits_{j=p+1}^{\infty}\int\limits_{\mathds{R}^{2}}\frac{e^{-\sqrt{j^{2}-p^{2}}|x_{2}-t_{2}|}}{\pi\sqrt{j^{2}-p^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}=0.$

In the same way we confirm that

$\displaystyle\operatorname{Im}\sum\limits_{j=1}^{p-1}\int\limits_{\mathds{R}^{2}}\frac{\tau\sin\tau\sqrt{p^{2}-j^{2}}|x_{2}-t_{2}|}{\pi\sqrt{p^{2}-j^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}=0,$
$\displaystyle\operatorname{Im}\frac{1}{\pi}\int\limits_{\mathds{R}^{2}}|x_{2}-t_{2}|U_{p}(x_{2})U_{p}(t_{2})\,dx_{2}\,dt_{2}=0.$

Hence, by two above identities and (3.7),

$\displaystyle\operatorname{Im}\int\limits_{\Omega}$ $\displaystyle
V_{1}(x)U_{\tau}(x)\sin
px_{1}\,dx=\operatorname{Im}\sum\limits_{j=1}^{p-1}\int\limits_{\mathds{R}^{2}}\frac{\mathrm{i}\tau\cos\tau\sqrt{p^{2}-j^{2}}|x_{2}-t_{2}|}{\pi\tau\sqrt{p^{2}-j^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}$
$\displaystyle=$
$\displaystyle\tau\operatorname{Re}\sum\limits_{j=1}^{p-1}\int\limits_{\mathds{R}^{2}}\frac{\cos\sqrt{p^{2}-j^{2}}|x_{2}-t_{2}|}{\pi\sqrt{p^{2}-j^{2}}}U_{j}(x_{2})U_{j}(t_{2})\,dt_{2}\,dx_{2}$
$\displaystyle=$
$\displaystyle\tau\sum\limits_{j=1}^{p-1}\int\limits_{\mathds{R}^{2}}\frac{\cos\sqrt{p^{2}-j^{2}}(x_{2}-t_{2})}{\pi\sqrt{p^{2}-j^{2}}}\big{(}W_{1,j}(x_{2})W_{1,j}(t_{2})-W_{2,j}(x_{2})W_{2,j}(t_{2})\big{)}\,dt_{2}\,dx_{2}.$

By straightforward calculations, for an arbitrary compactly supported function
$W(x_{2})$ we obtain:

$\displaystyle\int\limits_{\mathds{R}^{2}}\cos\sqrt{p^{2}-j^{2}}(x_{2}-t_{2})W(x_{2})W(t_{2})\,dt_{2}\,dx_{2}=$
$\displaystyle\left(\int\limits_{\mathds{R}}W(x_{2})\cos\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}$
$\displaystyle+\left(\int\limits_{\mathds{R}}W(x_{2})\sin\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}.$

Hence, by two latter identities and (3.9),

$\displaystyle\tau\operatorname{Im}\int\limits_{\Omega}V_{1}(x)U_{\tau}(x)\sin
px_{1}\,dx=\frac{1}{\pi}\sum\limits_{j=1}^{p-1}\frac{1}{\sqrt{p^{2}-j^{2}}}$
$\displaystyle\left(\left(\int\limits_{\mathds{R}}W_{1,j}(x_{2})\cos\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}\right.$
$\displaystyle\left.\hphantom{\Bigg{(}}-\left(\int\limits_{\mathds{R}}W_{2,j}(x_{2})\sin\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}\right).$

The latter formula implies that the sign of its left hand side is the same for
both $\tau\in\\{-1,+1\\}$ and we can make this sign being $-1$ by choosing
appropriately $W_{2}$ once we fix $W_{1}$ satisfying (3.10), namely, we can
satisfy the condition

$\displaystyle\sum\limits_{j=1}^{p-1}\frac{1}{\sqrt{p^{2}-j^{2}}}$
$\displaystyle\left(\left(\int\limits_{\mathds{R}}W_{1,j}(x_{2})\cos\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}\right.$
(3.11)
$\displaystyle\left.\hphantom{\Bigg{(}}-\left(\int\limits_{\mathds{R}}W_{2,j}(x_{2})\sin\sqrt{p^{2}-j^{2}}x_{2}\,dx_{2}\right)^{2}\right)<0$

For instance, this can be done by letting $W_{2}=\alpha\tilde{W}_{2}$ with a
sufficient large $\alpha$, where $\tilde{W}_{2}=\tilde{W}_{2}(x)$ is a real
odd compactly supported function such that

$\sum\limits_{j=1}^{p-1}\frac{1}{\sqrt{p^{2}-j^{2}}}\left(\int\limits_{\Omega}\tilde{W}_{2}(x)\sin
px_{1}\sin jx_{1}\sin\sqrt{p^{2}-j^{2}}x_{2}\,dx\right)^{2}>0.$

Once conditions (3.10), (3.11) hold, Theorem 2.5 states that both poles
$k_{\varepsilon,\tau}$, $\tau\in\\{-1,+1\\}$ correspond to the eigenvalues
located in the vicinity of the internal threshold $\Lambda_{p}$.

The above analytic results, namely, the discussed asymptotic expansions,
approximate well the true eigenvalues and resonances for sufficiently small
$\varepsilon$. In order to demonstrate how small $\varepsilon$ is to be
chosen, we make some numerical computations.

For numerics we use

$W_{1}(x)=-\sum_{j=1}^{3}a_{j}\sin jx_{1}\cos\frac{x_{2}}{2},\qquad
W_{2}(x)=\sum_{j=1}^{3}b_{j}\sin jx_{1}\sin x_{2},$

where $a_{j}$ and $b_{j}$ are real coefficients, and we additionally let
$W_{i}(x)\equiv 0$ as $|x_{2}|>\pi$, $i=1,2$. The corresponding eigenvalue
problem is approximated using a second-difference numerical scheme with
Dirichlet boundary conditions at $x_{1}=0$ and $x_{1}=\pi$ and a decay
condition at $x_{2}\to\pm\infty$. In order to achieve a numerically efficient
approximation of decay condition at $x_{2}\to\pm\infty$, a quasi-equidistant
grid [33] is used with a step size gradually increasing towards
$x_{2}\to\pm\infty$. For small $\varepsilon\lessapprox 0.2$, where the
localization of the eigenfunctions in $x_{2}$-direction is extremely weak, and
an adequate approximation of the decay condition as $x_{2}\to\pm\infty$ is
practically impossible, we use the Neumann condition in order to approximate
slowly decaying oscillating tails of the eigenfunctions:

$\partial_{x_{2}}\Psi(x_{1},\pm X_{0})=0,\quad\text{where}\quad X_{0}\gg 1.$

Figure 1: (a) Eigenvalue $\Lambda_{1}(\varepsilon)$ emerging from the bottom
of the spectrum computed using only the first term in expansion (3.4) [blue
curve] and two terms in (3.4) [green curve]. Red points connected by red lines
are obtained from direct numerical solution of the eigenvalue problem. (b) The
real and imaginary parts of the eigenfunction at $\varepsilon=0.1$ plotted as
a function of $x_{2}$ for $x_{1}=\pi/2$. Here $a_{1}=1$, $a_{3}=4$,
$b_{1}=0.5$ and all other coefficients are zero. Figure 2: The real and
imaginary parts of the eigenvalue $\Lambda_{2}(\varepsilon)$ emerging from an
internal threshold in the essential spectrum computed using only the first
term in expansion (3.5) [blue curve] and two terms in (3.5) [green curve]. Red
points connected by red lines are obtained from direct numerical evaluation of
the spectrum. (c) The real and imaginary parts of the eigenfunction at
$\varepsilon=0.3$ plotted as a function of $x_{2}$ for $x_{1}=\pi/4$. Two
insets show the plots of the squared amplitude of the eigenfunction
$|\Psi|^{2}=(\mathrm{Re\,}\Psi)^{2}+(\mathrm{Im\,}\Psi)^{2}$, using larger
domains in the horizontal axes and linear (left inset) and logarithmic (right
inset) scales in the vertical axes. (d) Full plot of the modulus of the
eigenfunction. In all panels $a_{1}=1$, $b_{2}=3$, and all other coefficients
$a_{i}$ and $b_{i}$ are zero. Only the eigenvalue with a positive imaginary
part is shown. There also exists a complex-conjugate eigenvalue with a
negative imaginary part and a $\mathcal{PT}$-conjugate eigenfunction.

In Fig. 1 we plot the dependencies $\Lambda_{1}(\varepsilon)$ obtained from
the asymptotic expansions, when only the leading term is taken into account
and both terms are used, for the particular set of parameters $a_{1}=1$,
$b_{1}=0.1$, $a_{2}=b_{2}=0$. For $\varepsilon\lessapprox 0.2$, the agreement
between the analytical predictions and numerical results is rather good, while
for large values of $\varepsilon$ it is only qualitative. In Fig. 2 we plot
the same data for two eigenvalues bifurcating from the internal threshold with
$p=2$; here only the eigenvalue with positive imaginary part is shown. The
following set of parameters is used: $a_{1}=1$, $b_{1}=0$, $a_{2}=0$, and
$b_{2}=3$. Again, the numerical results are in a good agreement with
asymptotic expansions for weak perturbations, and in a qualitative agreement
for stronger ones. Real part of the eigenvalue is found to decrease with the
increase of $\varepsilon$. Nevertheless, for all values of $\varepsilon$ shown
in Fig. 2 the real part of eigenvalue $\Lambda_{2}(\varepsilon)$ is larger
than the lower edge of the essential spectrum. Therefore, in the optical
context, the corresponding eigenfunction, whose three-dimensional modulus plot
is shown in Fig. 2(d), indeed represents a non-Hermitian generalization of a
bound state in the continuum.

####  3.3.2 Two-dimensional Bose-Einstein condensate with parabolic trapping

In the previous example a choice of $\omega$ was not really important and the
discussed results are of a more general nature. For instance, we can choose
the operator $\mathcal{H}^{\prime}$ being a quantum harmonic oscillator.
Namely, let $\omega=\mathds{R}$ and

$\mathcal{H}^{\prime}:=-\frac{d^{2}\
}{dx_{1}^{2}}+x_{1}^{2}\quad\text{on}\quad\mathds{R}.$

Then $\Omega=\mathds{R}^{2}$ and the operator $\mathcal{H}$ becomes

$\mathcal{H}=-\Delta+x_{1}^{2}\quad\text{on}\quad\mathds{R}^{2}.$

As a perturbation, we choose $\mathcal{L}(\varepsilon)=V_{1}$, where $V_{1}$
is a complex-valued compactly supported potential on $\mathds{R}^{2}$. Here

$\Lambda_{p}=2p+{1},\qquad\psi_{p}(x_{1})=\frac{e^{-x^{2}/2}}{\sqrt{2^{p}p!\sqrt{\pi}}}H_{p}\left(x_{1}\right),\qquad
H_{n}(t):=(-1)^{n}e^{t^{2}}\frac{d^{n}\ }{dt^{n}}e^{-t^{2}},$

i.e., $H_{p}$ are Hermite polynomials, and $p=0,1,\ldots$. Then all
calculations and results from previous example can be easily reproduced, just
in all formulae the functions $\frac{\sqrt{2}}{\sqrt{\pi}}\sin jx_{1}$ are to
be replaced by the functions $\psi_{j}(x_{1})$ introduced above.

The described situation corresponds to a two-dimensional cloud of Bose-
Einstein condensate with a parabolic confinement in $x_{1}$ direction.
Assuming that the interparticle interactions are negligible, that is, the
condensate is effectively linear, we can model its dynamics by the Schröginer-
like equation, which in the theory of Bose-Einstein condensates is known as
Gross-Pitaevskii equation [59]:

$\mathrm{i}\partial_{t}\Phi+\Delta\Phi-x_{1}^{2}\Phi-\varepsilon
V_{\varepsilon}(x)\Phi=0,$ (3.12)

where $\Phi(x_{1},x_{2},t)$ stands for the macroscopic wavefunction of the
condensate. Again, for stationary states in the form
$\Phi=e^{-\mathrm{i}\lambda t}\Psi$, where $\lambda$ has the meaning of the
chemical potential, the problem is reduced to the described eigenvalue
problem. Positive and negative imaginary parts of the perturbation correspond
to the injecting the particles from an external source and absorption of the
particles, respectively.

####  3.3.3 Three-dimensional circular waveguide and three-dimensional Bose-
Einstein condensate

Here we consider two examples of a three-dimensional waveguide and a three-
dimensional quantum oscillator, which extend the above examples of the planar
waveguide and the harmonic oscillator adduced above.

In first example we deal with a circular waveguide assuming that
$\omega=\\{x^{\prime}=(x_{1},x_{2}):\,|x^{\prime}|<1\\}$ and the operator
$\mathcal{H}^{\prime}$ is introduced as the Schrödinger operator with a
radially symmetric potential subject to the Dirichlet condition:

$\mathcal{H}^{\prime}=-\Delta_{x^{\prime}}+V_{0}\quad\text{in}\quad\omega,\qquad
V_{0}=V_{0}(|x^{\prime}|).$

Then the domain $\Omega$ is a straight cylinder along the axis $x_{3}$ with
the cross-section $\omega$ and

$\mathcal{H}=-\Delta+V_{0}(|x^{\prime}|)\quad\text{in}\quad\Omega$

subject to the homogeneous Dirichlet condition. The operator
$\mathcal{H}^{\prime}$ has a purely discrete spectrum and thanks to the
assumed radial symmetricity for the potential $V_{0}$, it possesses double
eigenvalues $\Lambda_{p}=\Lambda_{p+1}>\Lambda_{1}$ such that the associated
eigenfunctions, orthonormalized in $L_{2}(\omega)$, read as
$\psi_{p}(x^{\prime})=\Psi(|x^{\prime}|)\cos s\theta$,
$\psi_{p+1}(x^{\prime})=\Psi(|x^{\prime}|)\sin s\theta$, where $s$ is some
fixed integer number, $\Psi$ is some real function, and $\theta$ is a polar
angle associated with $x^{\prime}$. Such eigenvalues are degenerate internal
thresholds in the essential spectrum. We define then a
$\mathcal{PT}$-symmetric potential $V_{1}$ by formula (3.8) and in addition,
we assume that $W_{1}$ is even in $x_{1}$. Then it is easy to see that the
corresponding matrix $\mathrm{M}_{1}$ is diagonal:

$\mathrm{M}_{1}=\begin{pmatrix}\mu_{1}&0\\\
0&\mu_{2}\end{pmatrix},\qquad\mu_{1}:=-\frac{1}{2}\int\limits_{\Omega}W_{1}\psi_{p}^{2}\,dx,\qquad\mu_{2}:=-\frac{1}{2}\int\limits_{\Omega}W_{1}\psi_{p+1}^{2}\,dx.$

Then we assume that $\mu_{1}\neq\mu_{2}$ and $\mu_{1}<0$, $\mu_{2}<0$; these
conditions can be easily satisfied by choosing appropriately $W_{1}$. Since
for both $\mu_{i}$ we have $q_{i}=1$, for each corresponding pole we can apply
asymptotic expansion (2.44):

$k_{i,\tau}(\varepsilon)=\varepsilon\mu_{i}+\frac{\varepsilon^{2}}{2}\int\limits_{\Omega}V_{1}(x)U_{i,\tau}(x)\psi_{i-1+p}\,dx+O(\varepsilon^{3}),$

where $U_{i,\tau}:=\mathcal{G}_{1,\tau}V_{1}\psi_{i-1+p}$ are solutions to the
boundary value problem

$(-\Delta-\Lambda_{p})U_{i,\tau}=V_{1}\psi_{i-1+p}-\psi_{i-1+p}\int\limits_{\Omega}V_{1}(t^{\prime},x_{3})\psi_{p}^{2}(t^{\prime})\,dt^{\prime}{dx_{3}}\quad\text{in}\quad\Omega,\qquad
U_{i,\tau}=0\quad\text{on}\quad\partial\Omega,$

and can be found by a separation of variables similar to (3.6). Then we can
reproduce calculations from the first example and obtain that choosing
appropriately the function $W_{2}$, we can satisfy conditions (2.27), (2.34)
for all poles $k_{i,\tau}$, $i=1,2$, $\tau\in\\{-1,+1\\}$ and this means that
these poles correspond to eigenvalues. In other words, this means that for
appropriately chosen $W_{1}$ and $W_{2}$, we can generate four different
simple eigenvalues of the operator $\mathcal{H}_{\varepsilon}$ in the vicinity
of the double eigenvalue $\Lambda_{p}$ of the operator $\mathcal{H}^{\prime}$
serving as an internal threshold in the essential spectrum.

A similar situation can be realized also in other three-dimensional models.
For instance, we can let $\omega=\mathds{R}^{2}$ and
$\mathcal{H}^{\prime}=-\Delta_{x^{\prime}}+V_{0}(|x^{\prime}|)$, where
$V_{0}=V_{0}(t)$ is some function growing unboundedly at infinity. Such
operator $\mathcal{H}^{\prime}$ again can have double eigenvalues with the
eigenfunctions of form $\Psi(|x^{\prime}|)\cos s\theta$ and
$\Psi(|x^{\prime}|)\sin s\theta$. Linear combinations $\Psi(|x^{\prime}|)\cos
s\theta\pm\mathrm{i}\Psi(|x^{\prime}|)\sin s\theta$ correspond to vortex
states with integer $s$ being the vorticity or topological charge. Therefore,
in the context of Bose-Einstein condensates, this mechanism can be potentially
applied for generation of localized in all three spatial dimensions vortex
rings (see e.g. [19]).

We also observe that in both discussed three-dimensional examples the operator
$\mathcal{H}^{\prime}$ can have not only double eigenvalues, but also ones of
higher multiplicities $n$. And in such cases, it is possible to find a
$\mathcal{PT}$-symmetric potential $V_{1}$ generating $2n$ eigenvalues of the
operator $\mathcal{H}_{\varepsilon}$ in the vicinity of the considered
multiple eigenvalue.

##  4 Meromorphic continuation

In this section we prove Theorem 2.1 on the meromorphic continuation of the
resolvent of the operator $\mathcal{H}_{\varepsilon}$. First we prove some
auxiliary statements in a separate subsection and then we prove the theorem.

###  4.1 Auxiliary lemmata

In this subsection we prove three auxiliary lemmata, which will be employed in
the proof of Theorem 2.1. The first statement is Lemma 2.1.

###### Proof of Lemma 2.1.

Reproducing literally the proof of Lemma 2.3 in [9], one can check easily the
identity
$\operatorname{\sigma_{ess}}(\mathcal{H}_{\varepsilon})=\operatorname{\sigma_{ess}}(\mathcal{H})$.
The operator $\mathcal{H}$ can be represented as a sum of tensor products
$\mathcal{H}=\mathcal{H}^{\prime}\otimes\mathcal{I}+\mathcal{I}\otimes\mathcal{H}_{0}$.
And since

$\operatorname{\sigma_{ess}}(\mathcal{H}_{0})=\operatorname{\sigma}(\mathcal{H}_{0})=[0,+\infty),\qquad\inf\operatorname{\sigma}(\mathcal{H}^{\prime})=\Lambda_{1},$

we immediately get:

$\operatorname{\sigma}(\mathcal{H})=\operatorname{\sigma_{ess}}(\mathcal{H})=[\Lambda_{1},+\infty).$

This completes the proof. ∎

The next statement is Lemma 2.2.

###### Proof of Lemma 2.2.

For each $u\in\mathfrak{D}(\mathcal{H})\cap L^{\bot}$, each $\psi_{j}$ and
almost each $x_{d}\in\mathds{R}$ we have:

$\displaystyle\big{(}(\mathcal{H}^{\bot}u)(\cdot,x_{d}),\psi_{j}\big{)}_{L_{2}(\omega)}=$
$\displaystyle\mathfrak{h}^{\prime}\big{(}u(\cdot,x_{d}),\psi_{j}\big{)}-\left(\frac{\partial^{2}u}{\partial
x_{d}^{2}}(\cdot,x_{d}),\psi_{j}\right)_{L_{2}(\omega)}$ $\displaystyle=$
$\displaystyle-\Lambda_{j}\big{(}u(\cdot,x_{d}),\psi_{j}\big{)}_{L_{2}(\omega)}-\frac{d^{2}\
}{dx_{d}^{2}}\big{(}u(\cdot,x_{d}),\psi_{j}\big{)}_{L_{2}(\omega)}=0.$

Hence, the operator $\mathcal{H}^{\bot}$ maps $\mathfrak{D}(\mathcal{H})\cap
L^{\bot}$ into $L^{\bot}$. This is an unbounded self-adjoint operator in
$L^{\bot}$ associated with the restriction of the form $\mathfrak{h}$ on
$\mathfrak{D}(\mathfrak{h})\cap L^{\bot}$. Moreover, for each
$u\in\mathfrak{D}(\mathfrak{h})\cap L^{\bot}$ we have:

$\displaystyle\mathfrak{h}(u,u)=$
$\displaystyle\int\limits_{\mathds{R}}\mathfrak{h}^{\prime}\big{(}u(\cdot,x_{d}),u(\cdot,x_{d})\big{)}\,dx_{d}+\int\limits_{\Omega}\left|\frac{\partial
u}{\partial
x_{d}}\right|^{2}\,dx\geqslant\int\limits_{\mathds{R}}\mathfrak{h}^{\prime}\big{(}u(\cdot,x_{d}),u(\cdot,x_{d})\big{)}\,dx_{d}$
$\displaystyle\geqslant$ $\displaystyle
c_{0}\int\limits_{\mathds{R}^{d}}\|u(\cdot,x_{d})\|_{L_{2}(\omega)}^{2}\,dx_{d}=c_{0}\|u\|_{L_{2}(\Omega)}^{2}.$

Hence, the spectrum of the operator $\mathcal{H}^{\bot}$ is located in
$[c_{0},+\infty)$. The proof is complete. ∎

The third lemma provides a meromorphic continuation for the resolvent of the
unperturbed operator.

######  Lemma 4.1.

Fix $p\in\\{1,\ldots,m\\}$ and $\tau\in\\{-1,+1\\}$ and let
$\Lambda_{p}=\ldots=\Lambda_{p+n-1}$ be an $n$-multiple eigenvalue of the
operator $\mathcal{H}^{\prime}$, where $n\geqslant 1$. There exists a
sufficiently small fixed $\delta>0$ such that for all complex $k\in
B_{\delta}$ and all $f\in L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ the boundary
value problem

$\begin{gathered}\left(-\sum\limits_{i,j=1}^{d-1}\frac{\partial\ }{\partial
x_{i}}A_{ij}\frac{\partial\ }{\partial
x_{j}}+\mathrm{i}\sum\limits_{j=1}^{d-1}\left(A_{j}\frac{\partial\ }{\partial
x_{j}}+\frac{\partial\ }{\partial
x_{j}}A_{j}\right)+A_{0}-\Lambda_{p}+k^{2}\right)u=f\quad\text{in}\quad\Omega,\\\
\mathcal{B}u=0\quad\text{on}\quad\partial\Omega,\end{gathered}$ (4.1)

is solvable in $W_{2,loc}^{2}(\Omega)$ and possesses a solution, which can be
represented as

$u=\mathcal{A}_{1,\tau}(k)f,$ (4.2)

where $\mathcal{A}_{1,\tau}$ is a linear operator mapping
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ into
$W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$. This operator is bounded and
meromorphic in $k\in B_{\delta}$. It has the only pole in $B_{\delta}$, which
is at zero and simple:

$\displaystyle\mathcal{A}_{1,\tau}(k)=\frac{1}{k}\mathcal{A}_{2}+\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k),$
(4.3)
$\displaystyle\mathcal{A}_{2}f:=\sum\limits_{j=p}^{p+n-1}\psi_{j}\ell_{j}f,\qquad\ell_{j}f:=\frac{1}{2}\int\limits_{\Omega}\overline{\psi_{j}(x^{\prime})}f(x)\,dx,$
(4.4)

where $\mathcal{A}_{4,\tau}(k):\,L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)\to
W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ is an operator bounded uniformly
in $k\in\overline{B_{\delta}}$ and holomorphic in $k\in B_{\delta}$, and, we
recall, $\mathcal{G}_{p,\tau}$ is the operator defined in (2.10).

For sufficiently large $x_{d}$, the solution $u$ given by (4.2) can be
represented as

$u(x,k)=\sum\limits_{j=1}^{m}u_{j}^{\pm}(x_{d},k)\psi_{j}(x_{d})+u_{\bot}^{\pm}(x,k),\qquad\pm
x_{d}>R,$ (4.5)

where $R$ is some fixed number, $u_{j}^{\pm}\in
L_{2}(I_{\pm},e^{{\mp\vartheta}x_{d}}dx_{d})$ are some meromorphic in $k\in
B_{\delta}$ functions possessing the asymptotic behavior

$\displaystyle u_{j}^{\pm}(x_{d})=e^{-\tau
K_{j}(k)|x_{d}|}\big{(}C_{j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},\quad$
$\displaystyle|x_{d}|\to\infty,\qquad j=1,\ldots,p-1,$ (4.6) $\displaystyle
u_{j}^{\pm}(x_{d})=e^{-K_{j}(k)|x_{d}|}\big{(}C_{j}^{\pm}(k)+O({e}^{-{\tilde{\vartheta}}|x_{d}|})\big{)},$
$\displaystyle|x_{d}|\to\infty,\qquad j=p,\ldots,m,$

$C_{j}^{\pm}(k)$ are some meromorphic in $k\in B_{\delta}$ functions,
$0<{\tilde{\vartheta}}<{\vartheta}$ is some fixed constant independent of $k$
and $x$, and $u_{\bot}^{\pm}\in W_{2}^{2}(\Omega_{R}^{\pm})$ are some
functions meromorphic in $k\in B_{\delta}$ and obeying the identities

$(u_{\bot}^{\pm}(\cdot,x_{d}),\psi_{j})_{L_{2}(\Omega)}=0$ (4.7)

for almost each $x_{d}\in I_{\pm}$ and for each $j=1,\ldots,m$.

###### Proof.

The functions

$f_{j}(x_{d}):=\big{(}f(\cdot,x_{d}),\psi_{j}\big{)}_{L_{2}(\omega)}$

obviously belong to $L_{2}(\mathds{R},e^{{\vartheta}|x_{d}|}dx_{d})$ since
$f\in L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$. It is straightforward to check
that the functions

$\displaystyle u_{j}(x_{d},k):=\frac{1}{2\tau
K_{j}(k)}\int\limits_{\mathds{R}}e^{-\tau
K_{j}(k)|x_{d}-y_{d}|}f_{j}(y_{d})\,dy_{d}\quad$ $\displaystyle\text{as}\quad
j<p-1,$ (4.8) $\displaystyle
u_{j}(x_{d},k):=\frac{1}{2K_{j}(k)}\int\limits_{\mathds{R}}e^{-K_{j}(k)|x_{d}-y_{d}|}f_{j}(y_{d})\,dy_{d}\quad$
$\displaystyle\text{as}\quad j\geqslant p,$

solve the equations

$-u_{j}^{\prime\prime}+(\Lambda_{j}-\Lambda_{p}+k^{2})u_{j}=f_{j}\quad\text{in}\quad\mathds{R},\qquad
j=1,\ldots,m.$

Employing these facts, we seek a solution to problem (4.1) as

$u(x,k)=\sum\limits_{j=1}^{m}u_{j}(x_{d}^{\prime},k)\psi_{j}(x^{\prime})+u_{\bot}(x,k)$
(4.9)

and for $u_{\bot}$ we immediately get problem (4.1) with $f$ replaced by the
function $f^{\bot}$ defined in (2.11). We see easily that $f^{\bot}\in
L^{\bot}$.

Thanks to Lemma 2.2, the resolvent
$(\mathcal{H}^{\bot}-\Lambda_{p}+k^{2})^{-1}$ is well-defined for all
sufficiently small complex $k\in B_{\delta}$ provided $\delta$ is small
enough. This resolvent is holomorphic in $k\in B_{\delta}$ as an operator in
$L_{2}(\Omega)$ and is uniformly bounded in $k\in\overline{B_{\delta}}$. It
can be expanded via the standard Neumann series:

$(\mathcal{H}^{\bot}-\Lambda_{p}+k^{2})^{-1}=\sum\limits_{j=0}^{\infty}(-k^{2})^{j}\big{(}(\mathcal{H}^{\bot}-\Lambda_{p})^{-1}\big{)}^{j+1}.$
(4.10)

This expansion implies that the resolvent
$(\mathcal{H}^{\bot}-\Lambda_{p}+k^{2})^{-1}$ is also holomorphic as an
operator from $L_{2}(\Omega)$ into $W_{2}^{2}(\Omega)$. We define then

$u_{\bot}:=(\mathcal{H}^{\bot}-\Lambda_{p}+k^{2})^{-1}f^{\bot}$

and this obviously gives a solution to problem (4.1). We denote the operator
mapping $f$ into the described solution by $\mathcal{A}_{1,\tau}(k)$ and let
us show that it possesses all stated properties.

First of all we observe that the introduced operator is independent of the
choice of $\tau$ as $\Lambda_{p}=\Lambda_{1}$ since in this case the functions
$u_{j}$ with $\Lambda_{j}<\Lambda_{p}$ are missing and the above constructions
become independent of $\tau$.

Just by the embeddings $L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)\subset
L_{2}(\Omega)$ and $W_{2}^{2}(\Omega)\subset
W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$, the resolvent
$(\mathcal{H}^{\bot}-\Lambda_{p}+k^{2})^{-1}$ is a bounded operator from
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ into
$W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ holomorphic in $k\in B_{\delta}$
for sufficiently small $\delta$ and bounded uniformly in
$k\in\overline{B_{\delta}}$. An operator mapping $f$ into
$\sum\limits_{j=1}^{m}u_{j}(x_{d},k)\psi_{j}(x^{\prime})$ is given explicitly
and by straightforward calculations we can check that this is also a bounded
operator from $L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ into
$W_{2}^{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)$ holomorphic in $k\in
B_{\delta}$. The calculations are based on estimates of the following kind:

$\displaystyle\|u_{j}\psi_{j}\|_{L_{2}(\Omega,e^{-{\vartheta}|x_{d}|}dx)}^{2}=$
$\displaystyle\|u_{j}\|_{L_{2}(\mathds{R},e^{-{\vartheta}|x_{d}|}dx)}^{2}$
$\displaystyle\leqslant$
$\displaystyle\frac{\|f_{j}\|_{L_{2}(\mathds{R},e^{{\vartheta}|x_{d}|}dx_{d})}^{2}}{2|K_{j}|^{2}}\int\limits_{\mathds{R}^{2}}e^{-|x_{d}-y_{d}|\tau\operatorname{Re}K_{j}(k)-a(|x_{d}|+|y_{d}|)}\,dx_{d}\,dy_{d}$
$\displaystyle\leqslant$ $\displaystyle
C(k)\|f\|_{L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)}^{2},$

where $C(k)$ is some constant independent of $f$. Since the functions $f$
belong to $L_{2}(\Omega,e^{\vartheta|x_{d}|}dx_{d})$, the functions $u_{j}$
can be expanded into the Laurent series with respect to the small parameter
$k$; in particular, this means that the functionals
$\ell_{j}:L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)\to\mathds{C}$ are well-defined
and bounded. Substituting these expansions and (4.10) into (4.9), we
immediately get representation (4.3), (4.4), (2.10).

In view of formula (4.9), in order to prove representation (4.5), (4.6),
(4.7), it is sufficient to analyze the behavior of the functions $u_{j}$ at
infinity. This can be done easily by an obvious identity

$\displaystyle u_{j}(x_{d},k)-\frac{e^{\mp\tau K_{j}(k)x_{d}}}{2\tau
K_{j}(k)}\int\limits_{\mathds{R}}e^{\pm\tau
K_{j}(k)y_{d}}f_{j}(y_{d})\,dy_{d}=\tilde{u}_{j}^{\pm}(x_{d},k),\qquad\pm
x_{d}>0,$ (4.11)
$\displaystyle\tilde{u}_{j}^{\pm}(x_{d},k):=\frac{1}{2K_{j}(k)}\int\limits_{x_{d}}^{\pm\infty}\Big{(}e^{K_{j}(k)(x_{d}-y_{d})}-e^{-K_{j}(k)(x_{d}-y_{d})}\Big{)}f_{j}(y_{d})\,dy_{d},$

and an estimate

$\displaystyle|\tilde{u}_{j}^{\pm}(x_{d},k)|\leqslant$
$\displaystyle\frac{1}{2|K_{j}(k)|}\left(\int\limits_{x_{d}}^{\pm\infty}\Big{|}e^{K_{j}(k)(x_{d}-y_{d})}-e^{-K_{j}(k)(x_{d}-y_{d})}\Big{|}^{2}e^{-a|y_{d}|}\,dy_{d},\right)^{\frac{1}{2}}\|f_{j}\|_{L_{2}(\mathds{R},e^{{\vartheta}|x_{d}|}dx_{d})}$
(4.12) $\displaystyle\leqslant$ $\displaystyle
C(k)e^{-{\tilde{\vartheta}}|y_{d}|}\|f_{j}\|_{L_{2}(\mathds{R},e^{{\vartheta}|x_{d}|}dx_{d})},\qquad\pm
x_{d}>0,$

where $C(k)>0$ and $0<{\tilde{\vartheta}<\vartheta}$ are some constants
independent of $x_{d}$ and $f$. The proof is complete. ∎

###  4.2 Meromorphic continuation for the resolvent of the perturbed operator

This section is devoted to the proof of Theorem 2.1. We begin with rewriting
the operator equation

$(\mathcal{H}_{\varepsilon}-\Lambda_{p}+k^{2})u_{\varepsilon}=f,$

as a boundary value problem (2.3). Then we denote

$g_{\varepsilon}:=f-\varepsilon\mathcal{L}(\varepsilon)u_{\varepsilon}$

and rewrite (2.3) as problem (4.1) with the function $f$ replaced by
$g_{\varepsilon}$. According Lemma 4.1, such problem is solvable and there
exists a solution given by formula (4.2):

$u_{\varepsilon}=\mathcal{A}_{1,\tau}(k)g_{\varepsilon}.$ (4.13)

We substitute this formula into corresponding boundary value problem (2.3) and
this leads to an operator equation in the space
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$:

$g_{\varepsilon}+\varepsilon\mathcal{L}(\varepsilon)\mathcal{A}_{1,\tau}(k)g_{\varepsilon}=f.$

Then we substitute representation (4.3), (4.4), (2.10) into this equation:

$g_{\varepsilon}+\frac{\varepsilon}{k}\sum\limits_{j=p}^{p+n-1}\big{(}\ell_{j}g_{\varepsilon}\big{)}\mathcal{L}(\varepsilon)\psi_{j}+\varepsilon\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}g_{\varepsilon}=f.$
(4.14)

Thanks to the properties of the operators $\mathcal{G}_{p,\tau}$ and
$\mathcal{A}_{4,\tau}(k)$ described in Lemma 4.1, the operator
$\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}$
is bounded uniformly in $\varepsilon$ and $k$ as an operator in
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ and is holomorphic in $k\in
B_{\delta}$. Hence, for sufficiently small $\varepsilon$, the operator
$(\mathcal{I}+\varepsilon\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)})$
is boundedly invertible and the inverse is holomorphic in $k\in B_{\delta}$ as
an operator in $L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)$; hereinafter the symbol
$\mathcal{I}$ stands for the identity mapping. We denote

$\mathcal{A}_{5,\tau}^{\varepsilon}(k):=\big{(}\mathcal{I}+\varepsilon\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}\big{)}^{-1}$
(4.15)

and apply this operator to equation (4.14). This gives:

$g_{\varepsilon}+\frac{\varepsilon}{k}\sum\limits_{j=p}^{p+n-1}\big{(}\ell_{j}g_{\varepsilon}\big{)}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{j}=\mathcal{A}_{5,\tau}^{\varepsilon}(k)f.$
(4.16)

Our next step is to find $\ell_{j}g_{\varepsilon}$. Once we do this, we shall
be able to find $g_{\varepsilon}$ from equation (4.14) as

$g_{\varepsilon}=-\frac{\varepsilon}{k}\sum\limits_{j=p}^{p+n-1}\big{(}\ell_{j}g_{\varepsilon}\big{)}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{j}+\mathcal{A}_{5,\tau}^{\varepsilon}(k)f$
(4.17)

and recover the function $u_{\varepsilon}$ by formula (4.13).

We apply the functionals $\ell_{i}$ to equation (4.16) with all $i$ such that
$\Lambda_{i}=\Lambda_{p}$. As a result, we arrive to a system of linear
equations for the vector $\ell
g_{\varepsilon}:=(\ell_{j}g_{\varepsilon})_{j=p,\ldots,p+n-1}$; this is a
vector-column. The system reads as

$\left(\mathrm{E}+\frac{\varepsilon}{k}\mathrm{M}_{\varepsilon,\tau}(k)\right)\ell
g_{\varepsilon}=F_{\varepsilon}(k).$ (4.18)

Here $\mathrm{E}$ is the unit matrix and $\mathrm{M}_{\varepsilon,\tau}(k)$ is
a square matrix with entries
$M_{\varepsilon}^{ij}(k):=\ell_{i}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{j}$,
where $i$ counts the rows and $j$ does the columns in the matrix
$\mathrm{M}_{\varepsilon,\tau}(k)$, while the symbol $F_{\varepsilon}(k)$
denotes a vector column with coordinates
$\ell_{i}\mathcal{A}_{5,\tau}^{\varepsilon}(k)f$, $i=p,\ldots,p+n-1$. As it
was shown in the proof of Lemma 4.1, the functionals
$\ell_{i}:L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)\to\mathds{C}$ are bounded and
since the operator $\mathcal{A}_{5,\tau}^{\varepsilon}(k)$ is a bounded one in
$L_{2}(\Omega,e^{\vartheta|x_{d}|}dx)$, the introduced matrix
$\mathrm{M}_{\varepsilon,\tau}(k)$ and vector $F_{\varepsilon}(k)$ are well-
defined.

In view of the holomorphy of the operator
$\mathcal{A}_{5,\tau}^{\varepsilon}$, the entries of the matrix
$\mathrm{M}_{\varepsilon,\tau}(k)$ are holomorphic in $k$ and this implies
that the matrix $\mathrm{M}_{\varepsilon,\tau}(k)$ is holomorphic in $k\in
B_{\delta}$. Hence, the determinant of the matrix
$k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k)$ is a holomorphic in
$k$ function and therefore, the matrix
$(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))^{-1}$ is well-
defined as a meromorphic in $k$ matrix. This allows us to solve system (4.18):

$\ell
g_{\varepsilon}=k(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))^{-1}F_{\varepsilon}(k)=k\big{(}\ell_{j}^{\varepsilon}(k)f\big{)}_{j=p,\ldots,p+n-1},\qquad\ell_{j}g_{\varepsilon}=k\ell_{j}^{\varepsilon}(k)f,$
(4.19)

where $\ell_{j}^{\varepsilon}(k)$ are functionals on
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}dx)$ meromorphic in $k\in B_{\delta}$.
Substituting these formulae into (4.17), we can find $g_{\varepsilon}$ and
determine $u_{\varepsilon}$ by formula (4.13):

$\displaystyle u_{\varepsilon}=$
$\displaystyle\sum\limits_{j=p}^{p+n-1}\big{(}\ell_{j}^{\varepsilon}f\big{)}\Big{(}\psi_{j}-\varepsilon\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{j}\Big{)}+\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}\mathcal{A}_{5,\tau}^{\varepsilon}(k)f$
(4.20) $\displaystyle=$
$\displaystyle\sum\limits_{j=p}^{p+n-1}\big{(}\ell_{j}^{\varepsilon}f\big{)}\big{(}\mathcal{I}+\varepsilon\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}\mathcal{L}(\varepsilon)\big{)}^{-1}\psi_{j}+\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)\big{)}\mathcal{A}_{5,\tau}^{\varepsilon}(k)f=:\mathcal{R}_{\varepsilon,\tau}(k)f.$

The function $u_{\varepsilon}$ solves problem (2.3) and satisfies
representation (2.4), (2.5), (2.6). The latter are implied by similar
representations (4.5), (4.6), (4.7) and formula (4.13). Since for sufficiently
small $k$ the functions $K_{j}$, $\Lambda_{j}<\Lambda_{p}$, satisfy asymptotic
identities (2.43), the functions $e^{-\tau K_{j}(k)|x_{d}|}$,
$\Lambda_{j}<\Lambda_{p}$ and $e^{-k|x_{d}|}$ decay exponentially at infinity
as $\operatorname{Re}k>0$ and $\tau\operatorname{Im}k^{2}<0$. Hence, under the
same conditions, we have
$u_{\varepsilon}=(\mathcal{H}_{\varepsilon}-\Lambda_{p}+k^{2})^{-1}f$ and this
means that the operator $\mathcal{R}_{\varepsilon,\tau}$ provides a
meromorphic continuation of the resolvent
$(\mathcal{H}_{\varepsilon}-\Lambda_{p}+k^{2})^{-1}$.

We observe that the introduced operator $\mathcal{R}_{\varepsilon,\tau}$ is
meromorphic in $k\in B_{\delta}$ as an operator from
$L_{2}(\Omega,e^{{\vartheta}|x_{d}|}\,dx)$ into $W_{2}^{2}(\Omega)\cap
L_{2}(\Omega,e^{-{\vartheta}|x_{d}|}\,dx)$ due to formulae (4.13), (4.17),
(4.19) and representation (4.3), (4.4), (2.10).

According to formula (4.20), the poles of the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ defined in this formula coincide with the
poles of the functionals $\ell_{j}^{\varepsilon}(k)$. Let at least one of
these functionals have a pole at a point $k_{\varepsilon}\in B_{\delta}$. In
view of the definition of the functionals in (4.19), this means that the
matrix
$(k_{\varepsilon}\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k_{\varepsilon}))$
is degenerate and by the Cramer’s rule, system (4.18) with
$F_{\varepsilon}(k)=0$, $k=k_{\varepsilon}$ has a non-trivial solution, which
we denote by
$l^{\varepsilon}=(l^{\varepsilon}_{j})_{j:\,\Lambda_{j}=\Lambda_{p}}$. We
observe that by the definition,

$k_{\varepsilon}l_{i}^{\varepsilon}=-\varepsilon\sum\limits_{j=p}^{p+n-1}l_{j}^{\varepsilon}\ell_{i}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{j}.$
(4.21)

Employing these identities and (4.3), (4.4), it is straightforward to check
that the formula

$\displaystyle\psi_{\varepsilon}:=$
$\displaystyle-\varepsilon\lim\limits_{k\to
k_{\varepsilon}}\mathcal{A}_{1,\tau}(k)\sum\limits_{i=1}^{n}l_{i}^{\varepsilon}\mathcal{A}_{5,\tau}^{\varepsilon}(k)\mathcal{L}(\varepsilon)\psi_{i+p-1}$
(4.22) $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{n}l_{i}^{\varepsilon}\psi_{i+p-1}-\varepsilon\big{(}\mathcal{G}_{p,\tau}+k_{\varepsilon}\mathcal{A}_{4,\tau}(k_{\varepsilon})\big{)}\sum\limits_{i=1}^{n}l_{i}^{\varepsilon}\mathcal{A}_{5,\tau}^{\varepsilon}(k_{\varepsilon})\mathcal{L}(\varepsilon)\psi_{i+p-1}$
$\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{n}l_{i}^{\varepsilon}\big{(}\mathcal{I}+\varepsilon\big{(}\mathcal{G}_{p,\tau}+k_{\varepsilon}\mathcal{A}_{4,\tau}(k_{\varepsilon})\big{)}\mathcal{L}(\varepsilon)\big{)}^{-1}\psi_{i+p-1}$

defines a non-trivial solution to problem (2.3) with $f=0$,
$k=k_{\varepsilon}$. Thanks to formulae (4.9), (4.8), (4.12), this function
can be also represented as

$\displaystyle\psi_{\varepsilon}=\sum\limits_{j=1}^{m}\psi_{j}(x^{\prime})\varphi_{j}(x_{d})+(\mathcal{H}^{\bot}-\Lambda_{p}+k_{\varepsilon}^{2})^{-1}f^{\bot},$
(4.23) $\displaystyle\varphi_{j}(x_{d}):=\frac{1}{2\tau
K_{j}(k_{\varepsilon})}\int\limits_{\Omega}e^{-\tau
K_{j}(k_{\varepsilon})|x_{d}-y_{d}|}f(y)\overline{\psi_{j}(y^{\prime})}\,dy$
$\displaystyle\text{as}\quad j<p-1,$
$\displaystyle\varphi_{j}(x_{d}):=\frac{1}{2K_{j}(k_{\varepsilon})}\int\limits_{\Omega}e^{-K_{j}(k_{\varepsilon})|x_{d}-y_{d}|}f(y)\overline{\psi_{j}(y^{\prime})}\,dy$
$\displaystyle\text{as}\quad j\geqslant p+n-1,$
$\displaystyle\varphi_{j}(x_{d}):=\frac{1}{2k_{\varepsilon}}\int\limits_{\Omega}e^{-k_{\varepsilon}|x_{d}-y_{d}|}f(y)\overline{\psi_{j}(y^{\prime})}\,dy$
$\displaystyle\text{as}\quad j=p,\ldots,p+n-1,\quad k_{\varepsilon}\neq 0,$
$\displaystyle\varphi_{j}(x_{d}):=-\frac{1}{2}\int\limits_{\Omega}|x_{d}-y_{d}|f(y)\overline{\psi_{j}(y^{\prime})}\,dy$
$\displaystyle\text{as}\quad j=p,\ldots,p+n-1,\quad k_{\varepsilon}=0,$

with

$f:=-\varepsilon\sum\limits_{i=1}^{n}l_{i}^{\varepsilon}\mathcal{A}_{5,\tau}^{\varepsilon}(k_{\varepsilon})\mathcal{L}(\varepsilon)\psi_{i+p-1}$

and $f^{\bot}$ defined by (2.11). We also observe that as $k_{\varepsilon}=0$,
it follows from (4.21) that

$(f(\cdot,x_{d}),\psi_{j})_{L_{2}(\Omega)}=0$

for almost each $x_{d}\in I_{\pm}$ and for each $j=1,\ldots,m$.
Representations (4.22), (4.23) and relations (4.11), (4.12) imply
representation (2.7), (2.8), (2.9). This completes the proof of Theorem 2.1.

##  5 Poles of meromorphic continuations

In this section we study the asymptotic behaviour of poles of the operators
$\mathcal{R}_{\varepsilon,\tau}$ and we prove Theorems 2.2, 2.3. As it was
shown in the proof of Theorem 2.1 in the previous section, the poles of the
operator $\mathcal{R}_{\varepsilon,\tau}(k)$ coincides with poles of the
matrix $(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))^{-1}$. The
orders of the mentioned poles of the operator
$\mathcal{R}_{\varepsilon,\tau}(k)$ and the orders of the zeroes of
$\det(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))$ obviously
coincide.

###  5.1 Proof of Theorem 2.2

Let us study the solvability of the equation

$\det(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))=0$ (5.1)

in $B_{\delta}$. It is straightforward to check that

$\det(k\mathrm{E}+\varepsilon\mathrm{M}_{\varepsilon,\tau}(k))=k^{n}+P(k,\varepsilon),\qquad
P(k,\varepsilon):=\sum\limits_{i=1}^{n}\varepsilon^{i}k^{n-i}P_{i}(k,\varepsilon),$

where $P_{i}(k,\varepsilon)$ are some functions holomorphic in $k\in
B_{\delta}$ and uniformly bounded in $k\in\overline{B_{\delta}}$ and
sufficiently small $\varepsilon$. We fix arbitrary
$\delta^{\prime}\in(0,\delta)$ and for sufficiently small $\varepsilon$ we
have an uniform estimate:

$|P(k,\varepsilon)|\leqslant C\varepsilon\quad\text{on}\quad\partial
B_{\delta^{\prime}}.$ (5.2)

Since $|k^{n}|=(\delta^{\prime})^{n}$ on $\partial B_{\delta^{\prime}}$ and
the function $k\mapsto k^{n}$ has the only zero in $B_{\delta}$ of order $n$
at the origin, estimate (5.2) allows us to apply the Rouché theorem and to
conclude that equation (5.1) has exactly $n$ zeroes counting their orders and
these zeroes converge to the origin as $\varepsilon\to+0$.

Our next step is to show that all zeroes of equation (5.1) are located in a
ball $B_{b\varepsilon}$ for some fixed $b$. For arbitrary $b$ and sufficiently
small $\varepsilon$ such that $b\varepsilon<d$, by the definition of the
function $P$ we have:

$|P(k,\varepsilon)|\leqslant
C\varepsilon^{n}(1+b+\ldots+b^{n-1}),\qquad|k|^{n}=\varepsilon^{n}b^{n}\quad\text{on}\quad\partial
B_{b\varepsilon}.$

Then we choose $b$ large enough and $\varepsilon$ small enough so that

$b^{n}>C(1+b+\ldots+b^{n-1}),\qquad b\varepsilon<\delta,$

and applying the Rouché theorem once again, we see that the zeroes of equation
(5.1) are located inside the ball $B_{b\varepsilon}$. This fact allows us to
seek the zeroes of equation (5.1) as $k=z\varepsilon$, where $|z|<b$ for all
sufficiently small $\varepsilon$.

The operator $\mathcal{A}_{5,\tau}^{\varepsilon}$ defined in (4.15) can be
expanded into the standard Neumann series and this yields:

$\mathcal{A}_{5,\tau}^{\varepsilon}(k)=\mathcal{I}-\varepsilon\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k))+\varepsilon^{2}(\mathcal{L}(\varepsilon)\big{(}\mathcal{G}_{p,\tau}+k\mathcal{A}_{4,\tau}(k)))^{2}\mathcal{A}_{5,\tau}^{\varepsilon}(k).$
(5.3)

In view of the above identity, the change $k=z\varepsilon$, $|z|<b$, and the
definition of the matrix $\mathrm{M}_{{\varepsilon,\tau}}$, we can represent
the latter as

$\mathrm{M}_{\varepsilon,\tau}(k)=-\mathrm{M}_{1}+\varepsilon\tilde{\mathrm{M}}_{2}(z,\varepsilon),$
(5.4)

where $\tilde{\mathrm{M}}_{2}(z,\varepsilon)$ is some uniformly bounded matrix
holomorphic in $z$ and we recall that the matrix $\mathrm{M}_{1}$ is defined
in (2.12). Both these matrices are well-defined since they arise as the terms
in the expansion for $\mathrm{M}_{\varepsilon,\tau}$.

Then we can rewrite equation (5.1) as

$\det\big{(}z\mathrm{E}-\mathrm{M}_{1}+\varepsilon\tilde{\mathrm{M}}_{2}(z,\varepsilon)\big{)}=0.$
(5.5)

We recall that $\mu_{j}$, $j=1,\ldots,N$, are different eigenvalues of the
matrix $\mathrm{M}_{1}$ of the multiplicities $q_{j}$ and hence,

$\det(z\mathrm{E}-\mathrm{M}_{1})=\prod\limits_{j=1}^{N}(z-\mu_{j})^{q_{j}}.$

This allows us to represent equation (5.5) as

$\prod\limits_{i=1}^{N}(z-\mu_{i})^{q_{i}}+\varepsilon\tilde{P}(z,\varepsilon)=0,$
(5.6)

where $\tilde{P}$ is some function holomorphic in $z$ and bounded uniformly in
$z$ and $\varepsilon$:

$|\tilde{P}(z,\varepsilon)|\leqslant C.$ (5.7)

By $T_{r}(z_{0})$ we denote a ball of radius $r$ in the complex plane centered
at the point $z_{0}$. We introduce the balls
$T_{\tilde{b}\varepsilon^{\frac{1}{q_{i}}}}(\mu_{i})$, where $\tilde{b}$ is
some fixed number. Then for sufficiently small $\varepsilon$, each ball
$T_{\tilde{b}\varepsilon^{\frac{1}{q_{i}}}}(\mu_{i})$ contains the only zero
of the function $z\mapsto\prod\limits_{j=1}^{N}(z-\mu_{j})^{q_{j}}$, this zero
is $\mu_{i}$ and the order of this zero is $q_{i}$. On the boundary of each
ball $T_{\tilde{b}\varepsilon^{\frac{1}{q_{i}}}}(\mu_{i})$, for sufficiently
small $\varepsilon$, an obvious estimate holds true:

$\left|\prod\limits_{j=1}^{N}(z-\mu_{j})^{q_{j}}\right|\geqslant
C\tilde{b}^{q_{i}}\varepsilon,$ (5.8)

where $C$ is some fixed constant independent of $\varepsilon$ and $\tilde{b}$.
Hence, in view of estimate (5.7), we can choose $\tilde{b}$ large enough and
$\varepsilon$ small enough and apply the Rouché theorem to the left hand side
of equation (5.6). This yields that each ball
$T_{\tilde{b}\varepsilon^{\frac{1}{q_{i}}}}(\mu_{i})$ contains exactly $q_{i}$
zeroes of equation (5.6), which we denote by $z_{ij}=z_{ij}(\varepsilon)$.
These zeroes satisfy the asymptotic identities

$z_{ij}(\varepsilon)=\mu_{i}+O(\varepsilon^{\frac{1}{q_{i}}}),\qquad
j=1,\ldots,q_{i}.$ (5.9)

Returning back to equation (5.1), we see that it has exactly $n$ zeroes
counting their orders and these zeroes obey asymptotic expansion (2.13). The
proof is complete.

###  5.2 Proof of Theorem 2.3

Identity (5.3) allows us to specify the form of the matrix
$\mathrm{M}_{\varepsilon,\tau}$ in more details than in (5.4). Namely, this
identity implies that

$\mathrm{M}_{\varepsilon,\tau}(k)=-\mathrm{M}_{1}+\varepsilon\mathrm{M}_{2}^{\tau}+\varepsilon^{2}\tilde{\mathrm{M}}_{3}(z,\varepsilon),$

where $\tilde{\mathrm{M}}_{3}(z,\varepsilon)$ is some uniformly bounded matrix
holomorphic in $z$ and we recall that the matrix $\mathrm{M}_{2}$ is defined
in (2.14). All these matrices are well-defined as the terms in the expansion
for $\mathrm{M}_{\varepsilon,\tau}$. Equation (5.5) is replaced by a more
detailed one:

$\det\big{(}z\mathrm{E}-\mathrm{M}_{1}+\varepsilon\mathrm{M}_{2}^{\tau}+\varepsilon^{2}\tilde{\mathrm{M}}_{3}(z,\varepsilon)\big{)}=0.$
(5.10)

We fix $i\in\\{1,\ldots,N\\}$ and we are going to study the asymptotic
behavior of the zeroes $z_{ij}(\varepsilon)$, $j=1,\ldots,q_{i}$. Let
$\mathrm{S}$ be a matrix reducing the matrix $\mathrm{M}_{1}$ to its Jordan
canonical form, which we denote by
$\mathrm{J}:=\mathrm{S}\mathrm{M}_{1}\mathrm{S}^{-1}$. Then equation (5.10)
can be rewritten as

$\det\big{(}z\mathrm{E}-\mathrm{J}+\varepsilon\mathrm{S}\mathrm{M}_{2}^{\tau}\mathrm{S}^{-1}+\varepsilon^{2}\mathrm{S}\tilde{\mathrm{M}}_{3}(z,\varepsilon)\mathrm{S}^{-1}\big{)}=0.$

Taking into consideration the structure of the matrix $\mathrm{J}$, we rewrite
this equation as follows:

$\prod\limits_{j=1}^{N}(z-\mu_{j})^{q_{j}}+\varepsilon(z-\mu_{i})^{r_{i,\tau}}Y_{1}(z)+\varepsilon^{2}Y_{2}(z,\varepsilon)=0,$
(5.11)

where $Y_{1}$, $Y_{2}$ are some functions holomorphic in $z$ and bounded
uniformly in $\varepsilon$ and $z$ and $r_{i,\tau}<q_{i}$ is some non-negative
integer number. These functions arise as some polynomials in $z$,
$\varepsilon$ and the entries of the matrices
$\mathrm{S}\mathrm{M}_{2}\mathrm{S}^{-1}$ and
$\mathrm{S}\tilde{\mathrm{M}}_{3}(z,\varepsilon)\mathrm{S}^{-1}$. It is clear
that $(z-\mu_{i})^{r_{i,\tau}}Y_{1}(z)=Q_{i,\tau}(z)$, where $Q_{i}$ is
defined in (2.15). And if $Y_{1}$ is not identically zero, then

$Y_{1}(z_{i})=\gamma_{i,\tau}\prod\limits_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}(\mu_{i}-\mu_{j})^{q_{j}}\neq 0,$ (5.12)

where $\gamma_{i,\tau}$ is defined in (2.17). By $T_{r}(z_{0})$ we denote a
ball of radius $r$ in the complex plane centered at the point $z_{0}$.

We first consider the case, when $Y_{1}$ vanishes identically. Then the term
$\varepsilon Y_{1}$ disappears in equation (5.11). Here we consider the ball
$T_{\tilde{b}\varepsilon^{\frac{2}{q_{i}}}}(\mu_{i})$ and for sufficiently
small $\varepsilon$, on the boundary of this ball we have the estimates

$\Big{|}\prod\limits_{j=1}^{N}(z-\mu_{j})^{q_{j}}\Big{|}\geqslant
C\tilde{b}^{q_{i}}\varepsilon^{2},\qquad|\varepsilon^{2}Y_{2}(z,\varepsilon)|\leqslant\tilde{C}\varepsilon^{2},$

where $C$ and $\tilde{C}$ are some fixed constants independent of
$\varepsilon$ and $\tilde{b}$. Then we apply the Rouché theorem proceeding as
in (5.6), (5.7), (5.8), (5.9) and we arrive immediately to (2.16).

Now we proceed to the case, when $Y_{1}$ is not identically zero. Here we
divide equation (5.11) by $\prod\limits_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}(z-\mu_{j})^{q_{j}}$ and in view of (5.12) we get:

$(z-\mu_{i})^{q_{i}}+\varepsilon(z-\mu_{i})^{r_{i,\tau}}\gamma_{i,\tau}+\varepsilon(z-\mu_{i})^{r_{i,\tau}+1}Y_{3}(z)+\varepsilon^{2}Y_{4}(z,\varepsilon)=0,$
(5.13)

where $Y_{3}$, $Y_{3}$ are some functions holomorphic in $z$ in the vicinity
of the point $\mu_{i}$ and bounded uniformly in $\varepsilon$ and $z$.

Suppose that $q_{i}\leqslant 2r_{i,\tau}$ and consider the ball
$T_{\tilde{b}\varepsilon^{\frac{1}{r_{i,\tau}}}}(\mu_{i})$. For sufficiently
small $\varepsilon$, on the boundary of this ball we have

$\displaystyle|(z-\mu_{i})^{q_{i}}+\varepsilon(z-\mu_{i})^{r_{i,\tau}}\gamma_{i,\tau}|\geqslant$
$\displaystyle|(z-\mu_{i})^{q_{i}}|-\varepsilon|\gamma_{i,\tau}||(z-\mu_{i})^{r_{i,\tau}}|$
(5.14) $\displaystyle=$
$\displaystyle\tilde{b}^{q_{i}}\varepsilon^{2-\frac{2r_{i,\tau}-q_{i}}{r_{i,\tau}}}-\varepsilon^{2}|\gamma_{i,\tau}|\tilde{b}^{r_{i,\tau}}\geqslant\frac{\tilde{b}^{q_{i}}}{2}\varepsilon^{2-\frac{2r_{i,\tau}-q_{i}}{r_{i,\tau}}}$

and