  *

# Equations of motion governing the dynamics of the exceptional points of
parameterically dependent nonhermitian Hamiltonians

Milan Šindelka  Institute of Plasma Physics of the Czech Academy of Sciences,
Za Slovankou 1782/3, 18200 Prague, Czech Republic<EMAIL_ADDRESS>Pavel
Stránský and Pavel Cejnar  Institute of Nuclear and Particle Physics, Faculty
of Mathematics and Physics, Charles University, V Holešovičkách 2, 18000
Prague, Czech Republic<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study exceptional points (EPs) of a nonhermitian Hamiltonian
$\hat{H}(\lambda,\delta)$ whose parameters $\lambda\in{\mathbb{C}}$ and
$\delta\in{\mathbb{R}}$. As the real control parameter $\delta$ is varied, the
$k$-th EP (or $k$-th cluster of simultaneously existing EPs) of
$\hat{H}(\lambda,\delta)$ moves in the complex plane of $\lambda$ along a
continuous trajectory, $\lambda_{k}(\delta)$. We derive a self contained set
of equations of motion (EOM) for the trajectory $\lambda_{k}(\delta)$, while
interpreting $\delta$ as the propagation time. Such EOM become of interest
whenever one wishes to study the response of EPs to external perturbations or
continuous parametric changes of the pertinent Hamiltonian. This is e.g. the
case of EPs emanating from hermitian curve crossings/degeneracies (which turn
into avoided crossings/near-degeneracies when the Hamiltonian parameters are
continuously varied). The presented EOM for EPs have not only their
theoretical merits, they possess also a substantial practical relevance.
Namely, the just presented approach can be regarded even as an efficient
numerical method, useful for generating EPs for a broad class of complex
quantum systems encountered in atomic, nuclear and condensed matter physics.
Performance of such a method is tested here numerically on a simple yet
nontrivial toy model.

††: J. Phys. A: Math. Gen.

Keywords:
nonhermitian degeneracies, dynamics of exceptional points, avoided crossings.

##  1 Introduction

Nonhermitian Hamiltonians give rise to a special kind of degeneracies (the so
called exceptional points, EPs) which are not encountered within the standard
hermitian quantum mechanics. Namely, not only the (complex) eigenvalues, but
also the corresponding eigenvectors become degenerate (coalescent) at the EP
[1, 2, 3, 4, 5, 6]. Mathematical peculiarities of such a situation include
self-orthogonality, an unusual closure property, and multivaluedness of the
involved eigenvalues when encircling an EP in the parameter space of the
Hamiltonian. Importantly, the EPs arise not only in toy models, but also in a
vast amount of physically relevant and experimentally accessible contexts
(quantum mechanics of laser driven atoms, waveguide optics, acoustics,
electric circuit theory, elasticity) where they imply surprising counter-
intuitive phenomena, see e.g. short reviews [7, 8, 9] and also Refs. [10, 11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33]. Moreover, relevance of EPs in the context of quantum chaos and
quantum phase transitions has been demonstrated theoretically [11, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44]. The role of EPs in the superradiance
phenomenon has also been recognized [45, 46]. Higher order EPs have been
explored e.g. in Refs. [47, 48, 49].

Substantial effort has been invested into developing computational methods for
finding the EPs explicitly for a given nonhermitian Hamiltonian [50, 51, 52,
53]. In spite of great ingenuity and insighfulness of such algorithms, their
application to concrete systems is not always straightforward. Difficulties
arise especially when the Hamiltonian under study supports and existence of
many EPs, which are typically associated with avoided crossings encountered in
the framework of the pertinent hermitian theory.

The purpose of our present article is to further contribute both to the theory
and computational methodology related to EPs. Namely, we study the response of
EPs to continuous changes of the Hamiltonian parameters. Our intention is then
to calculate EPs of a given problem by means of a continuous parametric
propagation, starting from an arrangement where the EPs are trivial (or at
least easy) to find.

Our basic idea can be sketched as follows. Let us consider a parameterically
$\lambda$-dependent Hamiltonian of the general form

$\hat{H}(\lambda)\;=\;\hat{H}_{0}\;+\;\lambda\,\hat{V}\hskip 14.22636pt.\hskip
14.22636pt(\lambda\in{\mathbb{C}})$ (1)

We are looking for the EPs of $\hat{H}(\lambda)$ in the complex
$\lambda$-plane. The present paper pursues the following strategy: We
conveniently express $\hat{V}$ as a sum $\hat{V}=\hat{V}_{0}+\hat{V}_{1}$,
where the component $\hat{V}_{0}$ is chosen so that the eigenvalue problem and
the EPs of $\hat{H}_{0}+\lambda\hat{V}_{0}$ are either trivially resolvable or
at least easy to handle. A prototypical example (which will be elaborated
fully explicitly below in Section 3) corresponds to cases when
$[\hat{H}_{0},\hat{V}_{0}]=\hat{0}$. Then $\hat{H}_{0}$ and $\hat{V}_{0}$
possess the same eigenvectors, and situations closely linked to the EPs are
encountered due to the exact level crossings of
$\hat{H}_{0}+\lambda\hat{V}_{0}$. Our original Hamiltonian (1) can be now
redisplayed as

$\hat{H}(\lambda)\;=\;\hat{H}_{0}\;+\;\lambda\,\hat{V}_{0}\;+\;\lambda\,\hat{V}_{1}\hskip
14.22636pt.$ (2)

Formula (2) motivates us to think of a slightly more general Hamiltonian

$\hat{H}(\lambda,\delta)\;=\;\hat{H}_{0}\;+\;\lambda\,\hat{V}_{0}\;+\;\lambda\,\delta\,\hat{V}_{1}\hskip
14.22636pt;$ (3)

where $\delta\in[0,1]$ serves as an auxiliary switching parameter of the
$\lambda\,\hat{V}_{1}$ term. Importantly, one has
$\hat{H}(\lambda,0)=\hat{H}_{0}+\lambda\hat{V}_{0}$ and
$\hat{H}(\lambda,1)=\hat{H}(\lambda)$ of Eq. (2). Moreover, the EPs of
$\hat{H}(\lambda,\delta)$ of Eq. (3) move continuously in the complex
$\lambda$-plane when the real valued switching control parameter $\delta$ is
set to increase gradually from 0 to 1. If so, it seems natural to examine the
possibility of finding the equations of motion (EOM) governing the ”flux” or
”dynamical propagation” of the mentioned EPs along the ”time coordinate”
$\delta\in[0,1]$. One may even anticipate that an explicit solution of such
EOM (where the initial conditions at $\delta=0$ are provided by the presumably
known EPs of $\hat{H}(\lambda,0)=\hat{H}_{0}+\lambda\hat{V}_{0}$) would lead
to finding the desired EPs of $\hat{H}(\lambda)=\hat{H}(\lambda,1)$. It is the
purpose of our present article to adequately explore both theoretical and
practical merits of the just sketched approach.

The just presented idea corresponds essentially to implementing adequately and
systematically the nonhermitian perturbation theory in the presence of EPs,
where the perturbation is invoked by parameteric shift
$\delta\mapsto\delta+{\rm d}\delta$. Let us mention in this context that the
merits of nonhermitian perturbation theory in the presence of an EP have been
recently exploited e.g. in Ref. [16] (see the corresponding Supplementary
material). We also point out that our approach is intimately related to the
Dyson-Pechukas theory of level dynamics (see e.g. Ref. [43]), which however
has been pursued so far just within the framework of hermitian Hamiltonian
formalism.

The paper is organized as follows. Section 2 provides a self contained
systematic fully explicit theoretical derivation of the sought EOM for the
EPs. As such, Section 2 represents the most important hardcore material to be
communicated by the present article. The issue of choosing appropriate initial
conditions for the EOM is conveniently relegated to Appendix A. Section 3
describes a conceptually simple yet certainly nontrivial toy model, intended
to serve as a relatively strict test of our obtained EOM. Symmetry of this toy
model (resulting in simultaneous existence of multiple EPs at particular
values of $\lambda$) is discussed. Subsequently, we present in Section 3 an
outcome of the numerical solution of our EOM for the aforementioned toy model,
in order to highlight suitability of our EOM method for practical computations
of the EPs. Finally, Section 4 contains the concluding remarks.

##  2 Mathematical formulation

###  2.1 Preliminaries

Let us consider an $N$-by-$N$ complex symmetric111An extension of our
considerations to general non-symmetric Hamiltonians is also possible and
relatively straightforward. However, in the present article we prefer to deal
only with symmetric Hamiltonians for the sake of maximum simplicity.
Hamiltonian matrix $\hat{H}(\lambda,\delta)$ depending upon two parameters
$\lambda\in{\mathbb{C}}$ and $\delta\in{\mathbb{R}}$. The Hamiltonian
$\hat{H}(\lambda,\delta)$ is acting in the linear space ${\mathbb{C}}^{N}$ of
$N$-component ket (column) vectors

$|v)\;=\;\left(\matrix{v_{1}\cr v_{2}\cr\bm{\vdots}\cr v_{N}}\right)\hskip
14.22636pt.$ (4)

The associated bra (row) vectors are simply
$(v|=(v_{1}\;v_{2}\;\cdots\;v_{N})$. The adequate scalar product (the so
called $c$-product) is defined by prescription

$(v|v^{\prime})\;=\;\sum_{n=1}^{N}\,v_{n}\,v^{\prime}_{n}\;=\;(v^{\prime}|v)\hskip
14.22636pt.$ (5)

Recall that the self overlap $(v|v)$ is generally complex valued, and
$(v|v)=0$ does not imply $|v)$ equal to the zero vector $|\emptyset)$, see
Chapter 9 of Ref. [4] for details.

Consistently with our motivational considerations outlined in the
Introduction, we shall hereafter assume that there exists a function
$\lambda(\delta)$ such that, for each $\delta\in{\mathbb{R}}$, an eigenproblem
of $\hat{H}(\lambda(\delta),\delta)$ gives rise to $M$ distinct binary EPs
$[\,$$N\geq 2M$, each binary EP is formed via coalescence of two eigenvectors
of $\hat{H}(\lambda\to\lambda(\delta),\delta)$$\,]$.222 In the case when
$\hat{H}(\lambda(\delta),\delta)$ does not possess any kind of symmetry, we
expect $M=1$. On the other hand, symmetries of
$\hat{H}(\lambda(\delta),\delta)$ might imply $M>1$ (see Section 3 for an
example). We define for later convenience

$\hat{H}(\delta)\;\equiv\;\hat{H}(\lambda(\delta),\delta)\hskip 14.22636pt;$
(6)

and also

$\hat{V}(\delta)\;\equiv\;{\rm
d}_{\delta}\,\hat{H}(\delta)\;=\;\partial_{\lambda}\,\hat{H}(\lambda(\delta),\delta)\,\bm{\dot{}}{\lambda}(\delta)\;+\;\partial_{\delta}\,\hat{H}(\lambda(\delta),\delta)\hskip
14.22636pt;$ (7)

where $\bm{\dot{}}{\lambda}(\delta)={\rm d}_{\delta}\,\lambda(\delta)$, with
${\rm d}_{\bullet}=\frac{{\rm d}}{{\rm d}\bullet}$ and
$\partial_{\bullet}=\frac{\partial}{\partial\bullet}$.

The eigenproblem of our interest looks then as follows. In accordance with our
above made assumption, for each $\delta\in{\mathbb{R}}$ there exist $M$ binary
EPs of $\hat{H}(\delta)$, satisfying

$\hat{H}(\delta)\,|\tilde{c}_{m}^{\delta})\;=\;\tilde{E}_{m}^{\delta}\,|\tilde{c}_{m}^{\delta})\hskip
14.22636pt,\hskip 14.22636pt1\leq m\leq M\hskip 14.22636pt;$ (8)

with obvious notations.333 The upper tilde superscript indicates here entities
associated inherently with the EPs, whereas all the non-EP entities are
conveniently left without tilde. In this manner we distinguish e.g. between
$\tilde{E}_{1}^{\delta}$ of Eq. (8) and $E_{1}^{\delta}$ of Eq. (9).  Besides
these $M$ EPs, there exist also $(N-2M)$ ordinary non-degenerate non-EP
eigenvectors of $\hat{H}(\delta)$, satisfying

$\hat{H}(\delta)\,|c_{j}^{\delta})\;=\;E_{j}^{\delta}\,|c_{j}^{\delta})\hskip
14.22636pt,\hskip 14.22636pt1\leq j\leq(N-2M)\hskip 14.22636pt.$ (9)

Since the just listed ensemble of $(N-M)$ Hamiltonian eigenvectors
$|\tilde{c}_{m}^{\delta})$ and $|c_{j}^{\delta})$ does not form a complete
basis set of ${\mathbb{C}}^{N}$, one needs to include into the game also $M$
complementary basis vectors (see Section 9.2 of Ref. [4]), satisfying

$\left(\hat{H}(\delta)\,-\,\tilde{E}_{m}^{\delta}\,\hat{1}\right)\,|\tilde{b}_{m}^{\delta})\;=\;f_{m}^{\delta}\,|\tilde{c}_{m}^{\delta})\hskip
14.22636pt,\hskip 14.22636pt1\leq m\leq M\hskip 14.22636pt.$ (10)

Here $f_{m}^{\delta}$ are nonzero coefficients arising from imposing suitable
normalization conventions444 See equation (18) below and the accompanying
discussion. Subsection 2.2 and Appendix A describe an unambigous gauge fixing
of $f_{m}^{\delta}$ and all related matters. for $|\tilde{c}_{m}^{\delta})$
and $|\tilde{b}_{m}^{\delta})$, other notations are again self explanatory.

The corresponding orthonormality relations take the following explicit
appearance:

$\displaystyle(c_{j}^{\delta}|c_{j^{\prime}}^{\delta})$ $\displaystyle=$
$\displaystyle\delta_{jj^{\prime}}\hskip 14.22636pt;$ (11)
$\displaystyle(c_{j}^{\delta}|\tilde{c}_{m}^{\delta})$ $\displaystyle=$
$\displaystyle 0\hskip 14.22636pt;$ (12)
$\displaystyle(c_{j}^{\delta}|\tilde{b}_{m}^{\delta})$ $\displaystyle=$
$\displaystyle 0\hskip 14.22636pt;$ (13)
$\displaystyle(\tilde{c}_{m}^{\delta}|\tilde{c}_{m^{\prime}}^{\delta})$
$\displaystyle=$ $\displaystyle 0\hskip 14.22636pt;$ (14)
$\displaystyle(\tilde{c}_{m}^{\delta}|\tilde{b}_{m^{\prime}}^{\delta})$
$\displaystyle=$ $\displaystyle\delta_{mm^{\prime}}\hskip 14.22636pt;$ (15)
$\displaystyle(\tilde{b}_{m}^{\delta}|\tilde{b}_{m^{\prime}}^{\delta})$
$\displaystyle=$ $\displaystyle 0\hskip 14.22636pt.$ (16)

Relations (14), (16) show that the eigenvectors $|\tilde{c}_{m}^{\delta})$ and
their complements $|\tilde{b}_{m}^{\delta})$ are self orthogonal and
normalized via (15) and (10), as opposed to the eigenvectors
$|c_{j}^{\delta})$ which are unit normalizable through (11). The pertinent
closure property is built up accordingly (see again Section 9.2 of Ref. [4]),
we have

$\sum_{j}\,|c_{j}^{\delta})(c_{j}^{\delta}|\;+\;\sum_{m}\,|\tilde{c}_{m}^{\delta})(\tilde{b}_{m}^{\delta}|\,+\,|\tilde{b}_{m}^{\delta})(\tilde{c}_{m}^{\delta}|\;=\;\hat{1}\hskip
14.22636pt;$ (17)

where $\hat{1}$ stands for an $N$-by-$N$ unit matrix.

The normalization of self-orthogonal vectors $|\tilde{c}_{m}^{\delta})$ and
$|\tilde{b}_{m}^{\delta})$ is not unambiguously fixed by the formulas (10),
(15), (17). Indeed, these relations are invariant with respect to rescalings

$|\tilde{c}_{m}^{\delta,{\rm
new}})\;=\;g_{m}^{\delta}\;|\tilde{c}_{m}^{\delta})\hskip 7.11317pt,\hskip
7.11317pt|\tilde{b}_{m}^{\delta,{\rm
new}})\;=\;\left(g_{m}^{\delta}\right)^{-1}|\tilde{b}_{m}^{\delta})\hskip
7.11317pt,\hskip 7.11317ptf_{m}^{\delta,{\rm
new}}\;=\;\left(g_{m}^{\delta}\right)^{-2}f_{m}^{\delta}\hskip 7.11317pt;$
(18)

where $g_{m}^{\delta}$ stands for any nozero complex valued factor.

Our above outlined formulas (8)-(17) indicate that full solution of an
eigenvalue problem of $\hat{H}(\delta)$ is determined by seven fundamental
entities

$\lambda(\delta)\hskip 7.11317pt,\hskip 7.11317pt\tilde{E}_{m}^{\delta}\hskip
7.11317pt,\hskip 7.11317pt|\tilde{c}_{m}^{\delta})\hskip 7.11317pt,\hskip
7.11317pt|\tilde{b}_{m}^{\delta})\hskip 7.11317pt,\hskip
7.11317ptE_{j}^{\delta}\hskip 7.11317pt,\hskip 7.11317pt|c_{j}^{\delta})\hskip
7.11317pt,\hskip 7.11317ptf_{m}^{\delta}\hskip 14.22636pt.$ (19)

The just displayed entities (19) depend continuously upon the parameter
$\delta\in{\mathbb{R}}$. An infinitesimal shift of $\delta$ changes our
Hamiltonian (6) into

$\hat{H}(\delta+{\rm d}\delta)\;=\;\hat{H}(\delta)\;+\;\hat{V}(\delta)\,{\rm
d}\delta\hskip 14.22636pt.$ (20)

This invokes the corresponding infinitesimal changes in the eigensolutions
(19). The associated rates of change

$\bm{\dot{}}{\lambda}(\delta)\hskip 7.11317pt,\hskip
7.11317pt\bm{\dot{}}{\tilde{E}_{m}^{\delta}}\hskip 7.11317pt,\hskip
7.11317pt|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\hskip 7.11317pt,\hskip
7.11317pt|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\hskip 7.11317pt,\hskip
7.11317pt\bm{\dot{}}{E}_{j}^{\delta}\hskip 7.11317pt,\hskip
7.11317pt|\bm{\dot{}}{c}_{j}^{\delta})\hskip 7.11317pt,\hskip
7.11317pt\bm{\dot{}}{f}_{m}^{\delta}$ (21)

are obtainable by examining how do the eigensolutions (19) respond to the
Hamiltonian perturbation $\hat{V}(\delta)\,{\rm d}\delta$ in equation (20). An
explicit analytic elaboration of such a perturbation theory is by no means
conventional or trivial, since the considered eigenproblems of
$\hat{H}(\delta)$ and $\hat{H}(\delta+{\rm d}\delta)$ do support $M$ binary
EPs, as highlighted above in (8)-(17). Nevertheless, the just mentioned task
is feasible to perform, and results in explicit analytic prescriptions for the
”velocities” (21) determining the ”dynamics” or ”motion” of the seven
fundamental eigensolution entities (19) in the flux of ”time”
$\delta\in{\mathbb{R}}$. These ”equations of motion for the EPs” (or briefly
EOM) are worked out in a self contained manner in the next Subsection 2.2,
which actually represents the most important ”hard core” material to be
communicated by the present paper. See the resulting equations (24), (26),
(31), (43), (44), (45), (46) below. Furthermore, an additional Appendix A
describes in a self contained fashion the construction of adequate initial
conditions for these EOM, corresponding to a frequently encountered situation
when the sought EPs emanate from hermitian curve crossings/degeneracies of
$\hat{H}(\lambda,\delta)$.

###  2.2 Equations of motion for the exceptional points

Assume that the seven fundamental entities (19) are known for a given
$\delta\in{\mathbb{R}}$. Let us derive now in a self contained manner explicit
analytic formulas for the corresponding (presumably unknown) derivatives (21).
These need to be expressed solely in terms of the known quantities (19).

####  2.2.1 An equation of motion for $\bm{\dot{}}{\lambda}(\delta)$

bla

Take equation (8) for a given value of $m$ ($1\leq m\leq M$). Differentiate
both sides with respect to $\delta$, as to get

$\hat{V}(\delta)\,|\tilde{c}_{m}^{\delta})\;+\;\hat{H}(\delta)\,|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;\bm{\dot{}}{\tilde{E}_{m}^{\delta}}\,|\tilde{c}_{m}^{\delta})\;+\;\tilde{E}_{m}^{\delta}\,|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\hskip
14.22636pt.$ (22)

Substitute (7), multiply subsequently by $(\tilde{c}_{m}^{\delta}|$ from the
left, and exploit the self orthogonality property (14) at $m^{\prime}=m$ which
implies also

$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;0\hskip
14.22636pt.$ (23)

This yields a compelling formula

$\bm{\dot{}}{\lambda}(\delta)\;=\;-\,\frac{(\tilde{c}_{m}^{\delta}|\partial_{\delta}\,\hat{H}(\lambda(\delta),\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{c}_{m}^{\delta}|\partial_{\lambda}\,\hat{H}(\lambda(\delta),\delta)|\tilde{c}_{m}^{\delta})}\hskip
14.22636pt;$ (24)

which represents perhaps the most important result of the present paper.
Outcome (24) should be regarged as an equation of motion for
$\bm{\dot{}}{\lambda}(\delta)$. The r.h.s. of (24) must be independent upon
$m$, as long as our assumption of having $M$ binary EPs holds. The
$m$-independence of (24) serves as an useful check of internal consistency in
our numerical calculations of Section 3 of the main text. From now on,
$\bm{\dot{}}{\lambda}(\delta)$ will be regarded as explicitly known (and
presumably finite555It is beyond the scope of the present article to examine
if (or under which circumstances) the denominator of (24) can ever become
zero.), and the same applies also for the perturbation $\hat{V}(\delta)$ of
Eq. (7).

####  2.2.2 Equations of motion for $\bm{\dot{}}{\tilde{E}_{m}^{\delta}}$,
$\bm{\dot{}}{E}_{j}^{\delta}$, and $\bm{\dot{}}{f}_{m}^{\delta}$, plus other
accompanying elaborations

bla

Take equation (22) and multiply from the left by
$(\tilde{c}_{m^{\prime}}^{\delta}|$ where $m^{\prime}\neq m$. Exploit
subsequently (14). This yields an overlap element

$(\tilde{c}_{m^{\prime}}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\hskip
14.22636pt.\hskip 14.22636pt[\,m^{\prime}\neq m\,]$ (25)

The denominator of (25) is nonsingular as long as the considered $M$ binary
EPs are distinct.

Take again (22) and multiply from the left by $(\tilde{b}_{m}^{\delta}|$.
Exploit subsequently (10), (23), and also (15) for $m^{\prime}=m$. This yields
the as yet unknown energy derivative

$\bm{\dot{}}{\tilde{E}_{m}^{\delta}}\;=\;(\tilde{b}_{m}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})\hskip
14.22636pt.$ (26)

This is the sought equation of motion for
$\bm{\dot{}}{\tilde{E}_{m}^{\delta}}$.

Take again (22) and multiply from the left by
$(\tilde{b}_{m^{\prime}}^{\delta}|$ where $m^{\prime}\neq m$. Exploit
subsequently (10) and (15) together with (25). This yields an overlap element

$(\tilde{b}_{m^{\prime}}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\;+\;f_{m^{\prime}}^{\delta}\,\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}\hskip
14.22636pt.\hskip 14.22636pt[\,m^{\prime}\neq m\,]$ (27)

The denominators are nonsingular for the same reason as in (25). In passing we
note that

$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m^{\prime}}^{\delta}})\;=\;-\,(\tilde{b}_{m^{\prime}}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\hskip
14.22636pt;$ (28)

valid as an immediate consequence of (15).

Take again (22) and multiply from the left by $(c_{j}^{\delta}|$ where $1\leq
j\leq(N-2\,M)$. Exploit subsequently (9) and (12). This yields an overlap
element

$(c_{j}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;\frac{(c_{j}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-E_{j}^{\delta}}\hskip
14.22636pt.$ (29)

The denominator is nonsingular as long as the considered $m$-th EP eigenvalue
$\tilde{E}_{m}^{\delta}$ does not coincide with the non-EP eigenvalues
$E_{j}$.

Proceeding further, take equation (9) for a given value of $j$ ($1\leq j\leq
N-2\,M$). Differentiate both sides with respect to $\delta$, as to get

$\hat{V}(\delta)\,|c_{j}^{\delta})\;+\;\hat{H}(\delta)\,|\bm{\dot{}}{c}_{j}^{\delta})\;=\;\bm{\dot{}}{E}_{j}^{\delta}\,|c_{j}^{\delta})\;+\;E_{j}^{\delta}\,|\bm{\dot{}}{c}_{j}^{\delta})\hskip
14.22636pt.$ (30)

Multiply subsequently from the left by $(c_{j}^{\delta}|$, exploit then (9)
and (11) for $j^{\prime}=j$. This yields the as yet unknown energy derivative

$\bm{\dot{}}{E}_{j}^{\delta}\;=\;(c_{j}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})\hskip
14.22636pt.$ (31)

This is the sought equation of motion for $\bm{\dot{}}{E}_{j}^{\delta}$.

Take again (30) and multiply from the left by $(c_{j^{\prime}}^{\delta}|$
where $j^{\prime}\neq j$. Exploit then (9) and (11). This yields an overlap
element

$(c_{j^{\prime}}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})\;=\;\frac{(c_{j^{\prime}}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-E_{j^{\prime}}^{\delta}}\hskip
14.22636pt.\hskip 14.22636pt[\,j^{\prime}\neq j\,]$ (32)

The denominator is nonsingular as long as the non-EP eigenvalues are non-
degenerate.

Take again (30) and multiply from the left by $(\tilde{c}_{m}^{\delta}|$.
Exploit then (8) and (12). This yields an overlap element

$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})\;=\;\frac{(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-\tilde{E}_{m}^{\delta}}\hskip
14.22636pt.$ (33)

The denominator is nonsingular as long as $E_{j}$ does not coincide with the
EP eigenvalues $\tilde{E}_{m}^{\delta}$.

Take again (30) and multiply from the left by $(\tilde{b}_{m}^{\delta}|$.
Exploit then (10) and (13) together with (33). This yields an overlap element

$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})\;=\;\frac{(\tilde{b}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-\tilde{E}_{m}^{\delta}}\;+\;f_{m}^{\delta}\,\frac{(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{(E_{j}^{\delta}-\tilde{E}_{m}^{\delta})^{2}}\hskip
14.22636pt.$ (34)

Again, the denominators are nonsingular as long as $E_{j}$ does not coincide
with the EP eigenvalues $\tilde{E}_{m}^{\delta}$. In passing we note that

$(c_{j}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;=\;-\,(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})\hskip
14.22636pt;$ (35)

valid as an immediate consequence of (13).

Proceeding further, take equation (10) for a given value of $m$ ($1\leq m\leq
M$). Differentiate both sides with respect to $\delta$, as to get

$\hat{V}(\delta)\,|\tilde{b}_{m}^{\delta})\;+\;\hat{H}(\delta)\,|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;-\;\bm{\dot{}}{\tilde{E}_{m}^{\delta}}\,|\tilde{b}_{m}^{\delta})\;-\;\tilde{E}_{m}^{\delta}\,|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;=\;\bm{\dot{}}{f}_{m}^{\delta}\,|\tilde{c}_{m}^{\delta})\;+\;f_{m}^{\delta}\,|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\hskip
14.22636pt.$ (36)

Multiply subsequently from the left by $(\tilde{b}_{m}^{\delta}|$, exploit
then (10) and (15) for $m^{\prime}=m$, as well as (16) for $m^{\prime}=m$.
This yields another important relation

$(\tilde{b}_{m}^{\delta}|\hat{V}(\delta)\,|\tilde{b}_{m}^{\delta})\;=\;\bm{\dot{}}{f}_{m}^{\delta}\;+\;f_{m}^{\delta}\,(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;-\;f_{m}^{\delta}\,(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\hskip
14.22636pt.$ (37)

Take again (36) and multiply from the left by
$(\tilde{b}_{m^{\prime}}^{\delta}|$ where $m^{\prime}\neq m$. Exploit
subsequently (10) and (15), (16), together with (27) and (28). This yields an
overlap element

$\displaystyle(\tilde{b}_{m^{\prime}}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;=\hskip
284.52756pt[\,m^{\prime}\neq m\,]$ $\displaystyle=$
$\displaystyle\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\;+\;f_{m^{\prime}}^{\delta}\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}\;-\;f_{m}^{\delta}\;\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}\;-\;2\,f_{m}^{\delta}\,f_{m^{\prime}}^{\delta}\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{3}}\hskip
7.11317pt.$

Much like pointed out before, the denominator is nonsingular as long as the
considered $M$ binary EPs are distinct.

Take again (36) and multiply from the left by $(\tilde{c}_{m}^{\delta}|$.
Exploit subsequently (8), (14) for $m=m^{\prime}$, (15) for $m=m^{\prime}$,
and (23). This yields
$\bm{\dot{}}{\tilde{E}_{m}^{\delta}}=(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})$
as already known from (26).

Take again (36) and multiply from the left by
$(\tilde{c}_{m^{\prime}}^{\delta}|$ where $m^{\prime}\neq m$. Exploit
subsequently (8), (14), (15), together with (25). This yields an overlap
element
$(\tilde{c}_{m^{\prime}}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$ in the
form exactly equal to an immediate consequence of (27) and (28).

Take again (36) and multiply from the left by $(c_{j}^{\delta}|$ where $1\leq
j\leq(N-2\,M)$. Exploit subsequently (9), (12), (13), together with (29). This
yields an overlap element
$(c_{j}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$ in the form exactly
equal to an immediate consequence of (34) and (35).

To complete all our technical elaborations regarding the overlap elements, we
need to specify the as yet undetermined quantities
$(c_{j}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})$,
$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$,
$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$,
$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$. Recall that
$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$ is already
fixed by (23). Clearly, property (11) for $j^{\prime}=j$ implies immediately

$(c_{j}^{\delta}|\bm{\dot{}}{c}_{j}^{\delta})\;=\;0\hskip 14.22636pt.$ (39)

Similarly, property (16) for $m^{\prime}=m$ yields immediately

$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;=\;0\hskip
14.22636pt.$ (40)

An appropriate discussion of
$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$ and
$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$ is a bit more
intriguing. Property (15) for $m^{\prime}=m$ yields immediately

$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;-\,(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\hskip
14.22636pt;$ (41)

hence it is sufficient to determine just
$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$. Importantly,
the self-orthogonal vectors $|\tilde{c}_{m}^{\delta})$ and
$|\tilde{b}_{m}^{\delta})$, as well as the factor $f_{m}^{\delta}$, have been
introduced in the main text only modulo the rescaling transformation (18). It
is a trivial matter to verify that the as yet arbitrary rescaling coefficients
$g_{m}^{\delta}$ can be always chosen in such a particular manner as to
arrange for having

$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})\;=\;0\;=\;(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\hskip
14.22636pt.$ (42)

This is our suitably chosen gauge fixing convention for
$(\tilde{b}_{m}^{\delta}|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$ and
$(\tilde{c}_{m}^{\delta}|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$. Having imposed
(42), equation (37) simplifies into a finalized equation of motion for
$\bm{\dot{}}{f}_{m}^{\delta}$, namely,

$\bm{\dot{}}{f}_{m}^{\delta}\;=\;(\tilde{b}_{m}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})\hskip
14.22636pt.$ (43)

For the sake of completeness and clarity, let us also point out here that
$(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})=0$, this is
equivalent to (24).

####  2.2.3 Equations of motion for $|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$,
$|\bm{\dot{}}{c}_{j}^{\delta})$, $|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$

bla

The closure property (17) combined with (23), (25), (27), (29), (42) provides
immediately the desired equation of motion for
$|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$. One has

$\displaystyle|\bm{\dot{}}{\tilde{c}_{m}^{\delta}})$ $\displaystyle=$
$\displaystyle\sum_{j}\,|c_{j}^{\delta})\;\frac{(c_{j}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-E_{j}^{\delta}}$
(44) $\displaystyle+$ $\displaystyle\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\,\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\;+\;\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\,f_{m^{\prime}}^{\delta}\,\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}$
$\displaystyle+$ $\displaystyle\sum_{m^{\prime}\neq
m}\,|\tilde{b}_{m^{\prime}}^{\delta})\,\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\hskip
14.22636pt.$

Similarly, the closure property (17) combined with (32), (33), (34), (39)
provides immediately the desired equation of motion for
$|\bm{\dot{}}{c}_{j}^{\delta})$. One has

$\displaystyle|\bm{\dot{}}{c}_{j}^{\delta})$ $\displaystyle=$
$\displaystyle\sum_{j^{\prime}\neq
j}\,|c_{j^{\prime}}^{\delta})\,\frac{(c_{j^{\prime}}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-E_{j^{\prime}}^{\delta}}$
(45) $\displaystyle+$
$\displaystyle\sum_{m}\,|\tilde{c}_{m}^{\delta})\,\frac{(\tilde{b}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-\tilde{E}_{m}^{\delta}}\;+\;\sum_{m}\,|\tilde{c}_{m}^{\delta})\,f_{m}^{\delta}\,\frac{(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{(E_{j}^{\delta}-\tilde{E}_{m}^{\delta})^{2}}$
$\displaystyle+$
$\displaystyle\sum_{m}\,|\tilde{b}_{m}^{\delta})\,\frac{(\tilde{c}_{m}^{\delta}|\hat{V}(\delta)|c_{j}^{\delta})}{E_{j}^{\delta}-\tilde{E}_{m}^{\delta}}\hskip
14.22636pt.$

Finally, the closure property (17) combined with (27), (28), (34), (35),
(2.2.2), (40), (42) provides immediately the desired equation of motion for
$|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})$. One has

$\displaystyle|\bm{\dot{}}{\tilde{b}_{m}^{\delta}})\;=$ (46) $\displaystyle=$
$\displaystyle\sum_{j}\,|c_{j}^{\delta})\;\frac{(c_{j}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-E_{j}^{\delta}}\;-\;f_{m}^{\delta}\,\sum_{j}\,|c_{j}^{\delta})\;\frac{(c_{j}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-E_{j}^{\delta})^{2}}$
$\displaystyle+$ $\displaystyle\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\;\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}$
$\displaystyle+$ $\displaystyle\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\;f_{m^{\prime}}^{\delta}\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}\;-\;\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\;f_{m}^{\delta}\;\frac{(\tilde{b}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}$
$\displaystyle\hskip 147.95424pt\;-\;\sum_{m^{\prime}\neq
m}\,|\tilde{c}_{m^{\prime}}^{\delta})\;2\,f_{m}^{\delta}\,f_{m^{\prime}}^{\delta}\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{3}}$
$\displaystyle+$ $\displaystyle\sum_{m^{\prime}\neq
m}\,|\tilde{b}_{m^{\prime}}^{\delta})\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{b}_{m}^{\delta})}{\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta}}\;-\;f_{m}^{\delta}\,\sum_{m^{\prime}\neq
m}\,|\tilde{b}_{m^{\prime}}^{\delta})\;\frac{(\tilde{c}_{m^{\prime}}^{\delta}|\hat{V}(\delta)|\tilde{c}_{m}^{\delta})}{(\tilde{E}_{m}^{\delta}-\tilde{E}_{m^{\prime}}^{\delta})^{2}}\hskip
14.22636pt.$

In summary, we have in hand now a self contained collection of seven mutually
coupled equations of motion (24), (26), (31), (43), (44), (45), (46) for the
derivatives (21). These EOM determine the flux of the seven fundamental
entities (19) in the ”time” $\delta$, and can be propagated numerically along
$\delta\in{\mathbb{R}}$ once appropriate inital conditions are specified. As
already pointed out above, the issue of initial conditions is addressed in
Appendix A.

##  3 Test in a simple toy model

###  3.1 Introducing the toy model

The general mathematical formalism introduced above in Section 2 (and
supplemented by Appendix A) will be tested below on a conceptually simple yet
quite nontrivial toy model. We deliberately choose here a quantum system whose
state space is finite dimensional, and which therefore gives rise to a finite
number of EPs. Moreover, an inherent symmetry of our toy model allows
simultaneous existence of mutiple EPs ($M>1$), as anticipated already in
Section 2, and as explained in detail in the Figures below.

Our considered toy model corresponds to a theory of two distinct mutually
coupled angular momenta, $\hat{\vec{I}}=(\hat{I}_{1},\hat{I}_{2},\hat{I}_{3})$
and $\hat{\vec{J}}=(\hat{J}_{1},\hat{J}_{2},\hat{J}_{3})$, which possess the
conventional commutation properties. The associated starting Hamiltonian is
defined through prescription

$\hat{H}(\lambda)\;=\;\hat{H}_{0}\;+\;\lambda\,\hat{V}\hskip 14.22636pt;$ (47)

with

$\hat{H}_{0}\;=\;\omega\,\Bigl{(}\hat{I}_{3}+\hat{J}_{3}\Bigr{)}\hskip
14.22636pt;$ (48)

and

$\hat{V}\;=\;\hat{I}_{+}\,\hat{J}_{-}\,+\,\hat{I}_{-}\,\hat{J}_{+}\,+\,\hat{I}_{+}\,\hat{J}_{+}\,+\,\hat{I}_{-}\,\hat{J}_{-}\;=\;4\,\hat{I}_{1}\,\hat{J}_{1}\hskip
14.22636pt.$ (49)

Here $\omega>0$ and $\lambda\in{\mathbb{C}}$, and of course

$\hat{I}_{\pm}\;=\;\hat{I}_{1}\pm i\,\hat{I}_{2}\hskip 14.22636pt,\hskip
14.22636pt\hat{J}_{\pm}\;=\;\hat{J}_{1}\pm i\,\hat{J}_{2}\hskip 14.22636pt.$
(50)

The pertinent state space is spanned by basis vectors $|\,I_{\rm
T}\,I_{3}\,J_{\rm T}\,J_{3}\,\rangle$, where $(I_{\rm T}(I_{\rm T}+1),I_{3})$
are eigenvalues of $(\hat{I}^{2},\hat{I}_{3})$ and similarly $(J_{\rm
T}(J_{\rm T}+1),J_{3})$ are eigenvalues of $(\hat{J}^{2},\hat{J}_{3})$.
Clearly, both $I_{\rm T}$ and $J_{\rm T}$ are good quantum numbers for the
Hamiltonian (47). Dimension of a particular $(I_{\rm T},J_{\rm T})$ sector
equals to ${\cal N}_{I}\,{\cal N}_{J}$, here ${\cal N}_{I}=2\,I_{\rm T}+1$,
and similarly for ${\cal N}_{J}$. Hereafter we shall assume for definiteness
$I_{\rm T}=\frac{N}{2}$ and $J_{\rm T}=\frac{1}{2}$, where $N$ is an odd
positive integer (correspondingly, ${\cal N}_{I}=N+1$ and ${\cal N}_{J}=2$).
The parity of $(I_{3}+J_{3})$ is then another good quantum number.

Our primary interest consists in finding all the EPs of $\hat{H}(\lambda)$ of
Eq. (47) in the complex $\lambda$-plane. To accomplish this goal, we follow
the general strategy outlined in the Introduction, include into the game an
auxiliary switching parameter $\delta\in[0,1]$, and focus on investigating the
parameterically $\delta$-dependent EPs of an augmented Hamiltonian (3), where
by definition

$\hat{V}_{0}\;=\;\hat{I}_{+}\,\hat{J}_{-}\,+\,\hat{I}_{-}\,\hat{J}_{+}\hskip
14.22636pt;$ (51)
$\hat{V}_{1}\;=\;\hat{I}_{+}\,\hat{J}_{+}\,+\,\hat{I}_{-}\,\hat{J}_{-}\hskip
14.22636pt.$ (52)

Written explicitly, we have

$\hat{H}(\lambda,\delta)\;=\;\omega\,\Bigl{(}\hat{I}_{3}+\hat{J}_{3}\Bigr{)}\;+\;\lambda\,\left\\{\hat{I}_{+}\,\hat{J}_{-}\,+\,\hat{I}_{-}\,\hat{J}_{+}\,+\,\delta\,\Bigl{(}\hat{I}_{+}\,\hat{J}_{+}\,+\,\hat{I}_{-}\,\hat{J}_{-}\Bigr{)}\right\\}\hskip
7.11317pt.$ (53)

Again, both $I_{\rm T}$ and $J_{\rm T}$ are good quantum numbers for the
Hamiltonian (53), and the parity of $(I_{3}+J_{3})$ is another good quantum
number.

Before proceeding further, let us highlight an additional symmetry of the
Hamiltonian $\hat{H}(\lambda,\delta)$ of Eq. (53). The eigenvalue spectrum of
$\hat{H}(\lambda,\delta)$ is clearly invariant with respect to any similarity
transformation. One such particular transformation is represented by an
unitary operator

$\hat{U}\;=\;e^{-i\pi\hat{I}_{1}}\;e^{-i\pi\hat{J}_{2}}\hskip 14.22636pt;$
(54)

which corresponds to rotation of $\hat{\vec{I}}$ by angle $\pi$ around the
first coordinate axis, and to rotation of $\hat{\vec{J}}$ by angle $\pi$
around the second coordinate axis. Direct calculation yields

$\hat{U}^{\dagger}\,\hat{I}_{1}\,\hat{U}\;=\;+\,\hat{I}_{1}\hskip
14.22636pt,\hskip
14.22636pt\hat{U}^{\dagger}\,\hat{I}_{2}\,\hat{U}\;=\;-\,\hat{I}_{2}\hskip
14.22636pt,\hskip
14.22636pt\hat{U}^{\dagger}\,\hat{I}_{3}\,\hat{U}\;=\;-\,\hat{I}_{3}\hskip
14.22636pt;$ (55)

and similarly

$\hat{U}^{\dagger}\,\hat{J}_{1}\,\hat{U}\;=\;-\,\hat{J}_{1}\hskip
14.22636pt,\hskip
14.22636pt\hat{U}^{\dagger}\,\hat{J}_{2}\,\hat{U}\;=\;+\,\hat{J}_{2}\hskip
14.22636pt,\hskip
14.22636pt\hat{U}^{\dagger}\,\hat{J}_{3}\,\hat{U}\;=\;-\,\hat{J}_{3}\hskip
14.22636pt.$ (56)

Hence

$\hat{U}^{\dagger}\,\hat{I}_{\pm}\,\hat{U}\;=\;+\,\hat{I}_{\mp}\hskip
14.22636pt,\hskip
14.22636pt\hat{U}^{\dagger}\,\hat{J}_{\pm}\,\hat{U}\;=\;-\,\hat{J}_{\mp}\hskip
14.22636pt.$ (57)

If so, then our Hamiltonian $\hat{H}(\lambda,\delta)$ of Eq. (53) is converted
into

$\hat{U}^{\dagger}\,\hat{H}(\lambda,\delta)\,\hat{U}\;=\;-\,\hat{H}(\lambda,\delta)\hskip
14.22636pt.$ (58)

The just derived symmetry property (58) reveals that both
$\hat{H}(\lambda,\delta)$ and $-\hat{H}(\lambda,\delta)$ must possess the same
spectrum. Thus, if $E(\lambda,\delta)$ is an eigenvalue, then also
$-E(\lambda,\delta)$ is an eigenvalue. Note that another unitary
transformation $\hat{U}=e^{-i\pi\hat{I}_{2}}\,e^{-i\pi\hat{J}_{1}}$ leads to
the same conclusion, since the Hamiltonian (53) is invariant under an
interchange $\hat{\vec{I}}\leftrightarrow\hat{\vec{J}}$.

The matrix elements of $\hat{H}(\lambda,\delta)$ in a given $(I_{\rm T},J_{\rm
T})$ sector can be trivially calculated with the aid of familiar formulas for
the involved angular momentum operators. Let us write them down here
explicitly for the sake of maximum clarity:

$\displaystyle\langle I_{\rm T}\,I_{3}\,J_{\rm
T}\,J_{3}|\,\hat{H}(\lambda,\delta)\,|I_{\rm T}\,I^{\prime}_{3}\,J_{\rm
T}\,J^{\prime}_{3}\rangle\;=$ $\displaystyle=$
$\displaystyle\delta_{I_{3}I^{\prime}_{3}}\,\delta_{J_{3}J^{\prime}_{3}}\;\omega\,\Bigl{(}I_{3}+J_{3}\Bigr{)}$
$\displaystyle+$
$\displaystyle\delta_{I_{3}(I^{\prime}_{3}+1)}\,\delta_{J_{3}(J^{\prime}_{3}-1)}\;\lambda\,\phantom{\delta}\,\gamma_{+}(I_{\rm
T},I^{\prime}_{3})\,\gamma_{-}(J_{\rm
T},J^{\prime}_{3})\;+\;\delta_{I_{3}(I^{\prime}_{3}-1)}\,\delta_{J_{3}(J^{\prime}_{3}+1)}\;\lambda\,\phantom{\delta}\,\gamma_{-}(I_{\rm
T},I^{\prime}_{3})\,\gamma_{+}(J_{\rm T},J^{\prime}_{3})$ $\displaystyle+$
$\displaystyle\delta_{I_{3}(I^{\prime}_{3}+1)}\,\delta_{J_{3}(J^{\prime}_{3}+1)}\;\lambda\,\delta\,\gamma_{+}(I_{\rm
T},I^{\prime}_{3})\,\gamma_{+}(J_{\rm
T},J^{\prime}_{3})\;+\;\delta_{I_{3}(I^{\prime}_{3}-1)}\,\delta_{J_{3}(J^{\prime}_{3}-1)}\;\lambda\,\delta\,\gamma_{-}(I_{\rm
T},I^{\prime}_{3})\,\gamma_{-}(J_{\rm T},J^{\prime}_{3})\hskip 14.22636pt;$

where by definition

$\gamma_{\pm}(l,l^{\prime})\;=\;\sqrt{l(l+1)-l^{\prime}(l^{\prime}\pm
1)}\hskip 14.22636pt.$ (60)

In the case of $\delta=0$, the sum $K=(I_{3}+J_{3})$ becomes another good
quantum number. Correspondingly, the $(I_{\rm T},J_{\rm T})$ sector is divided
into subsectors associated with $K=(-I_{\rm T}-J_{\rm T}),(-I_{\rm T}-J_{\rm
T}+1),\cdots,(+I_{\rm T}+J_{\rm T})$. Moreover, $\hat{H}_{0}$ commutes both
with $\hat{I}_{+}\,\hat{J}_{-}$ and with $\hat{I}_{-}\,\hat{J}_{+}$, hence

$\Bigl{[}\hat{H}_{0}\,,\,\hat{V}_{0}\Bigr{]}\;=\;\hat{0}\hskip 14.22636pt;$
(61)

exactly as mentioned in the Introduction. Thereby an eigenvalue problem of the
Hamiltonian

$\hat{H}(\lambda,0)\;=\;\hat{H}_{0}\;+\;\lambda\,\hat{V}_{0}$ (62)

is solvable trivially, provided only that an eigenproblem of $\hat{V}_{0}$ has
been resolved. Accordingly, all the EPs of $\hat{H}(\lambda,0)$ are trivially
known (see the $\delta=0$ panels of Figs. 1 and 2 below, which consist just of
intersecting straight lines). This confirms that our definition of the
augmented Hamiltonian (53) satisfies the general requirements imposed on
$\hat{H}(\lambda,0)$ in the Introduction and in Appendix A.

For illustration, let us present now explicitly the calculated eigenvalue
spectrum of $\hat{H}(\lambda,\delta)$ for $N=19$, $\lambda\in[0,1]$,
$\delta\in[0,1]$, even parity of $(I_{3}+J_{3})$, and $\omega=1.0$. The
obtained results are shown in Fig. 1. An analogous case of odd parity is then
depicted in Fig. 2.

Figure 1:  The calculated eigenvalue spectrum of $\hat{H}(\lambda,\delta)$ of
Eq. (53) for $N=19$, even parity of $(I_{3}+J_{3})$, and $\omega=1.0$.
Horizontal axis corresponds to $\lambda$, vertical axis to the energy variable
$E$ associated with the eigenvalues. Note the reflection symmetry of the
spectrum with respect to the horizontal $E=0$ axis. This kind of symmetry is
explained by equation (58) above. Figure 2:  The calculated eigenvalue
spectrum of $\hat{H}(\lambda,\delta)$ of Eq. (53) for $N=19$, odd parity of
$(I_{3}+J_{3})$, and $\omega=1.0$. Horizontal axis corresponds to $\lambda$,
vertical axis to the energy variable $E$ associated with the eigenvalues. Note
the reflection symmetry of the spectrum with respect to the horizontal $E=0$
axis. This kind of symmetry is explained by equation (58) above.

The sought EPs of our starting Hamiltonian $\hat{H}(\lambda)$ of Eq. (47) can
be identified now with the EPs of $\hat{H}(\lambda,\delta)$ of Eq. (53) at
$\delta=1$. Yet the EPs of $\hat{H}(\lambda,\delta)$ are obtainable
numerically from the hermitian straight line crossings of $\hat{H}(\lambda,0)$
of Eq. (62) via the parametric $\delta$-propagation
$(\delta=0\mapsto\delta=1)$ of the EOM, exactly as we formulated in a self
contained fashion in the above Section 2 and in Appendix A.

###  3.2 Numerical propagation of the EOM and the obtained results

The seven mutually coupled equations of motion (24), (26), (31), (43), (44),
(45), (46) derived in Subsection 2.2 are propagated numerically using the
simplest possible first order difference scheme, starting from the initial
conditions which are established in Appendix A. At each propagation step, an
internal consistency of the obtained results is strictly checked. Namely, the
seven entities (19) calculated for a given particular value of
$\delta\in[0,1]$ are required to satisfy (up to a prescribed numerical
accuracy) the three eigenvalue equations (8), (9), (10), the six
orthonormality relations (11)-(16), and the closure property (17). In this
manner our numerical results presented below are granted to be reliably
converged.

Our illustrative numerical calculations are performed for the toy model of
Subsection 3.1, assuming $N=19$ and $\omega=1$ much as in Figs. 1-2. The range
$[0,1]$ of $\delta$ is discretized by $G=10^{7}$ equidistant grid points. This
ensures that our aforementioned test relations (8), (9), (10), (11)-(16), (17)
are fulfilled at each value of $\delta$ with the maximum error not exceeding
$0.0005$.

####  3.2.1 Results for the odd parity

bla

For maximum clarity of the presentation, it is convenient to start with
discussing our results obtained for the case of odd parity. Our propagation
starts from the hermitian crossings which are indicated by red bullets in Fig.
3.

Figure 3:  The hermitian straight line crossings corresponding to
$\hat{H}(\lambda,0)$, again for $N=19$, odd parity of $(I_{3}+J_{3})$, and
$\omega=1.0$. One may observe that there often (through not always) exist
multiple (twofold) crossings for a given value of $\lambda$ (see the vertical
black dashed lines). These are exactly the multiplets described theoretically
in Appendix A.

Our explicit numerical propagation of the EOM provides the following outcomes:

  * $\star$

Each isolated (onefold) hermitian crossing of Fig. 3 provides for $\delta>0$
an isolated (onefold) EP.

  * $\star$

Each twofold hermitian crossing of Fig. 3 provides for $\delta>0$ the
corresponding pair (twofold cluster) of distinct binary EPs which share the
same dependence $\lambda(\delta)$. This is a direct consequence of the
symmetry of $\hat{H}(\lambda,\delta)$ which is highlighted by equation (58)
above.

Figs. 4, 5 and 6 present explicitly our most important numerical results,
namely, the trajectories of the EPs in the $\lambda$-plane and in the plane of
complex energy. Note that Fig. 4 and Fig. 5 display essentially the same data,
just with a different layout convention. Specifically, Fig. 4 is plotted using
a random (machine generated) sign convention for the imaginary part of each
obtained curve $\lambda_{k}(\delta)$, whereas Fig. 5 corresponds to imposing a
fixed convention of $\Im\lambda_{k}(\delta)\geq 0$.666 Recall in this context
that each curve $\lambda_{k}(\delta)$ displayed in Figs. 4-5 gives rise to
another legitimate curve $\lambda^{*}(\delta)$ which penetrates into the
opposite side of the imaginary $\lambda$-plane, this curve
$\lambda_{k}^{*}(\delta)$ is not plotted.  Compared to Fig. 4, the overall
appearance of Fig. 5 is somewhat less transparent. In particular, some curves
$\lambda_{k}(\delta)$ do intersect (albeit at mutually distinct values of
$\delta$). This is the only reason why we hereafter prefer to present all our
numerical results using the layout convention analogous to Fig. 4.

Figure 4:  The EP trajectories $\lambda_{k}(\delta)$ emanating from the
hermitian straight line crossings of Fig. 3. The dark blue trajectories
correspond to a pair (twofold cluster) of distinct binary EPs which share the
same $\lambda_{k}(\delta)$, see our discussion in the main text. On the other
hand, the light blue trajectories are associated with a single binary EP. Note
also that each curve $\lambda_{k}(\delta)$ gives rise to another legitimate
curve $\lambda_{k}^{*}(\delta)$, which departs from the same (cluster of) red
bullet(s) of Fig. 3, but which corresponds to the complex conjugated initial
conditions at $\delta=0$. $[\,$This means that each curve
$\lambda_{k}(\delta)$ plotted explicitly here in the present figure has been
obtained via adopting a particular (machine generated) sign convention for the
$\sigma_{1}$-factor from Appendix A.$\,]$ Figure 5:  The same data as in Fig.
4, just a sign convention of $\Im\lambda_{k}(\delta)\geq 0$ is imposed a
posteriori. Note that each curve $\lambda_{k}(\delta)$ displayed here gives
rise to another legitimate curve $\lambda_{k}^{*}(\delta)$ which penetrates
into the negative imaginary plane of $\lambda$. Compared to Fig. 4, the
overall appearance of the present figure is somewhat less transparent. In
particular, some curves $\lambda_{k}(\delta)$ do intersect (albeit at mutually
distinct values of $\delta$). This is the only reason why we hereafter prefer
to display all our numerical results using the layout convention analogous to
Fig. 4.  Figure 6:  The EP trajectories $\tilde{E}_{m}(\delta)$ emanating from
the hermitian straight line crossings of Fig. 3 and corresponding to all the
curves $\lambda_{k}(\delta)$ plotted explicitly in Fig. 4. Importantly, all
the onefold hermitian crossings of Fig. 3 are associated with
$\tilde{E}_{m}(0)=0$, and actually provide $\tilde{E}_{m}(\delta)=0$ for all
$\delta\in[0,1]$. This fact (arising as a trivial consequence of the symmetry
property (58) of $\hat{H}(\lambda,\delta)$) is highlighted by the presence of
red bullet at the origin of the energy plane. On the other hand, the present
figure depicts also a progression of several nonzero trajectories
$\tilde{E}_{m}(\delta)$, which possess reflection symmetry with respect to the
origin. Each pair of these symmetry related trajectories corresponds
inevitably to a pair (twofold cluster) of distinct binary EPs which share the
same $\lambda_{k}(\delta)$.

####  3.2.2 Results for the even parity

bla

Let us move now on to the case of even parity. Our propagation starts from the
hermitian crossings which are indicated by red bullets in Fig. 7.

Figure 7:  The hermitian straight line crossings corresponding to
$\hat{H}(\lambda,0)$, again for $N=19$, even parity of $(I_{3}+J_{3})$, and
$\omega=1.0$. One may observe that there often (through not always) exist
multiple (twofold, fourfold) crossings for a given value of $\lambda$ (see the
vertical black dashed lines). These are exactly the multiplets described
theoretically in Appendix A.

Our explicit numerical propagation of the EOM provides the following outcomes:

  * $\star$

Each isolated (onefold) hermitian crossing of Fig. 7 provides for $\delta>0$
an isolated (onefold) EP.

  * $\star$

Each twofold hermitian crossing of Fig. 7 provides for $\delta>0$ the
corresponding pair (twofold cluster) of distinct binary EPs which share the
same dependence $\lambda(\delta)$. This is a direct consequence of the
symmetry of $\hat{H}(\lambda,\delta)$ which is highlighted by equation (58)
above.

  * $\star$

The fourfold crossings behave for $\delta>0$ as two separate twofold
crossings. Each of these twofold crossings reflects again the symmetry
property (58) of $\hat{H}(\lambda,\delta)$.

Figs. 8 and 9 present explicitly our most important numerical results, namely,
the trajectories of the EPs in the $\lambda$-plane and in the plane of complex
energy. We use the same layout convention as in Fig. 4 above.

Figure 8:  The EP trajectories $\lambda_{k}(\delta)$ emanating from the
hermitian straight line crossings of Fig. 7. The dark blue trajectories
correspond to a pair (twofold cluster) of distinct binary EPs which share the
same $\lambda_{k}(\delta)$, see our discussion in the main text. On the other
hand, the light blue trajectories are associated with a single binary EP. Note
also that each curve $\lambda_{k}(\delta)$ gives rise to another legitimate
curve $\lambda_{k}^{*}(\delta)$, which departs from the same (cluster of) red
bullet(s) of Fig. 7, but which corresponds to the complex conjugated initial
conditions at $\delta=0$. $[\,$This means that each curve
$\lambda_{k}(\delta)$ plotted explicitly here in the present figure has been
obtained via adopting a particular (machine generated) sign convention for the
$\sigma_{1}$-factor from Appendix A.$\,]$ Figure 9:  The EP trajectories
$\tilde{E}_{m}(\delta)$ emanating from the hermitian straight line crossings
of Fig. 7 and corresponding to all the curves $\lambda_{k}(\delta)$ plotted
explicitly in Fig. 8. Importantly, all the onefold hermitian crossings of Fig.
7 are associated with $\tilde{E}_{m}(0)=0$, and actually provide
$\tilde{E}_{m}(\delta)=0$ for all $\delta\in[0,1]$. This fact (arising as a
trivial consequence of the symmetry property (58) of
$\hat{H}(\lambda,\delta)$) is highlighted by the presence of red bullet at the
origin of the energy plane. On the other hand, the present figure depicts also
a progression of several nonzero trajectories $\tilde{E}_{m}(\delta)$, which
possess reflection symmetry with respect to the origin. Each pair of these
symmetry related trajectories corresponds inevitably to a pair (twofold
cluster) of distinct binary EPs which share the same $\lambda_{k}(\delta)$.

Summarizing the contents of Section 3, we have employed a nontrivial toy model
to explicitly test the performance of our computational algorithm based upon
solving the EOM for the EPs. We hope that our illustrative calculations
demonstrate practical usefulness of our EOM method for finding the EPs of
nontrivial Hamiltonians.

##  4 Concluding remarks

In summary, the present article establishes the equations of motion (EOM)
governing the dynamics (or flux) of EPs of parameterically dependent
nonhermitian Hamiltonians. This motion of EPs in the parameter space is
triggered here by a continuous change of an additional external control
parameter of the Hamiltonian. Our analysis covers a relatively broad class of
problems (1), where the search for EPs can be reinterpreted as solution of EOM
pertaining to an augmented Hamiltonian $\hat{H}(\lambda,\delta)$ of Eq. (3),
with $\delta$ playing the role of the dynamical ”time”.

From the theoretical point of view, Section 2 represents the most important
new hardcore material brought in by our paper. The resulting EOM (24), (26),
(31), (43), (44), (45), (46) are based essentially upon implementing a
nontraditional perturbation theory of nonhermitian quantum mechanics in the
presence of multiple EPs. An elaboration of such EOM, and in particular
derivation of equation (24), brings further theoretical insights into the
properties of the EPs, and represents thus a contribution of its own right.

Furthermore, our EOM can be exploited even in a totally pragmatic fashion,
merely as an efficient numerical tool for obtaining all the EPs of interest
for a given Hamiltonian $\hat{H}(\lambda)$ of Eq. (1). Such an approach lends
itself for its immediate application e.g. whenever the sought EPs emanate from
avoided crossings of the particular hermitian Hamiltonian under study. Section
3 demonstrates very explicitly practical merits of our EOM method in the just
mentioned situation.

We hope that the EOM formalism developed here can motivate or facilitate
further studies of EPs in atomic, nuclear, optical and condensed matter
physics.

Acknowledgements

We acknowledge financial support of the Czech Science Foundation under grant
Nos. 20-21179S (M. Š.) and 20-09998S (P. S. and P. C.), and of the Charles
University in Prague under project UNCE/SCI/013 (P. S. and P. C.).

References

## References

  * [1] T. Kato, Perturbation Theory of Linear Operators, Springer, New York (1966).
  * [2] C. M. Bender and T. T. Wu, Phys. Rev., 184, 1231 (1969).
  * [3] N. Moiseyev and S. Friedland, Phys. Rev. A, 22, 618 (1980).
  * [4] N. Moiseyev, Non-Hermitian Quantum Mechanics, Cambridge University Press (2011).
  * [5] A. P. Seyranian, O. N. Kirillov, and A. Mailybaev, J. Phys. A: Math. Gen., 38, 1723 (2005).
  * [6] M. A. Miri and A. Al, Science, 363, eaar7709 (2019).
  * [7] W. D. Heiss, J. Phys. A: Math. Gen., 37, 2455 (2004).
  * [8] I. Rotter, J. Phys. A: Math. Theor., 42, 153001 (2009).
  * [9] W. D. Heiss, J. Phys. A: Math. Theor., 45, 444016 (2012).
  * [10] M. R. Zirnbauer, J. J. M. Verbaarschot, and H. A. Weidenmüller, Nucl. Phys. A, 411, 161 (1983).
  * [11] W. D. Heiss and A. L. Sannino, J. Phys. A: Math. Gen., 23, 1167 (1990).
  * [12] W. D. Heiss and W. H. Steeb, J. Math. Phys., 32, 3003 (1991).
  * [13] M. V. Berry, Czech. J. Phys., 54, 1039 (2004).
  * [14] S. Garmon, M. Gianfreda, and N. Hatano, Phys. Rev. A, 92, 022125 (2015).
  * [15] S. Klaiman, U. Günther, and N. Moiseyev, Phys. Rev. Lett., 101, 080402 (2008).
  * [16] A. Pick, P. R. Kaprálová-Žďánská, N. Moiseyev, J. Chem. Phys., 150, 204111 (2019).
  * [17] C. Shi, M. Dubois, et al., Nature Communications, 7, 11110 (2016).
  * [18] Y. Choi, C. Hahn, et al., Nature Communications, 9, 2182 (2018).
  * [19] R. El-Ganainy, K. G. Makris, et al., Nature Physics, 14, 11 (2018).
  * [20] G. Shmuel, N. Moiseyev, Phys. Rev. Applied, 13, 024074 (2020).
  * [21] A. Ben-Asher, D. Šimsa, T. Uhlířová, M. Šindelka, and N. Moiseyev,
Phys. Rev. Lett., 124, 253202 (2020).

  * [22] N. Moiseyev and M. Šindelka, Phys. Rev. A, 103, 033518 (2021).
  * [23] P. R. Kaprálová, M. Šindelka, and N. Moiseyev, J. Phys. A: Math. Gen., 55, 284001 (2022).
  * [24] P. R. Kaprálová, Annals of Physics, 443, 168939 (2022).
  * [25] H. Hodaei, M. A. Miri, M. Heinrich, D. N. Christodoulides, M. Khajavikhan,
Science, 346, 975978 (2014).

  * [26] L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, X. Zhang, Science, 346, 972975 (2014).
  * [27] A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, U. Peschel,
Nature, 488, 167 (2012).

  * [28] Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, D. N. Christodoulides,
Phys. Rev. Lett., 106, 213901 (2011).

  * [29] X. Wang, X. Fang, D. Mao, Y. Jing, Y. Li, Phys. Rev. Lett., 123, 214302 (2019).
  * [30] W. Chen, S. K. Özdemir, G. Zhao, J. Wiersig, L. Yang, Nature, 548, 192 (2017).
  * [31] M. P. Hokmabadi, A. Schumer, D. N. Christodoulides, M. Khajavikhan, Nature, 576, 70 (2019).
  * [32] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides,
M. Khajavikhan, Nature, 548, 187 (2017).

  * [33] J. Wiersig, Phys. Rev. Lett., 112, 203901 (2014).
  * [34] W. D. Heiss, Z. Phys. A, 329, 133 (1989).
  * [35] W. D. Heiss and M. Müller, Phys. Rev. E, 66, 016217 (2002).
  * [36] W. D. Heiss, F. G. Scholtz, and H. B. Geyer, J. Phys. A: Math. Gen., 38, 1843 (2005).
  * [37] P. Cejnar, S. Heinze, and J. Dobeš, Phys. Rev. C, 71, 011304 (2005).
  * [38] P. Cejnar, S. Heinze, and M. Macek, Phys. Rev. Lett., 99, 100601 (2007).
  * [39] P. Stránský, M. Dvořák, and P. Cejnar, Phys. Rev. E, 97, 012112 (2018).
  * [40] T. E. Lee, F. Reiter, and N. Moiseyev, Phys. Rev. Lett., 113, 250401 (2014).
  * [41] D. I. Borisov, F. Ružička, and M. Znojil, Int. J. Theor. Phys., 54, 4293 (2015).
  * [42] M. Znojil, Proc. R. Soc. A, 476, 20190831 (2020).
  * [43] S.-J. Wang and S. Y. Chu, Phys. Rev. A, 47, 3546 (1993).
  * [44] M. Šindelka, L. F. Santos, N. Moiseyev, Phys. Rev. A, 95, 010103(R) (2017).
  * [45] C. Jung, M. Müller, and I. Rotter, Phys. Rev. E, 60, 114 (1999).
  * [46] P. Stránský and P. Cejnar, Phys. Rev. E, 100, 042119 (2019).
  * [47] E.-M. Graefe, U. Günther, H. J. Korsch, and A. E. Niederle,
J. Phys. A: Math. Theor., 41, 255206 (2008).

  * [48] G. Demange and E.-M. Graefe, J. Phys. A: Math. Theor., 45, 025303 (2012).
  * [49] M. Znojil, Phys. Rev. A, 100, 032124 (2019).
  * [50] A. Mailybaev, Numer. Linear Algebra Appl., 13, 419 (2006).
  * [51] R. Uzdin and R. Lefebvre, J. Phys. B: At. Mol. Opt. Phys., 43, 235004 (2010).
  * [52] O. N. Kirillov, Entropy, 20, 502 (2018).
  * [53] B. Nennig, E. Perrey-Debain, J. Comp. Phys., 412, 109425 (2020).

##  Appendix A Initial conditions for the EOM

Equations of motion (24), (26), (31), (43), (44), (45), (46) need to be
supplemented with appropriate initial conditions (ICS), i.e., by the seven
fundamental entities (19) provided at some starting value of
$\delta=\delta_{\rm in}$. The choice of $\delta_{\rm in}$ is of course
governed by concrete nature of the problem under study. In the present paper,
we shall describe a relatively frequently encountered situation when the
mentioned ICS are determinable at $\delta_{\rm in}$ (semi)trivially due to a
particularly simple form of $\hat{H}(\lambda,\delta_{\rm in})$. Namely, we
shall be concerned with such an arrangement when the starting Hamiltonian
$\hat{H}(\lambda,\delta_{\rm in})$ ($\lambda\in{\mathbb{R}}$) of the studied
physical model is hermitian (actually, even real symmetric), and possesses
exact crossings (accidental degeneracies). These crossings play the role of
origins from which our sought EPs emanate into the complex $\lambda$-plane as
$\delta$ is set to depart continuously from $\delta_{\rm in}$.

Let us assume for now $\lambda\in{\mathbb{R}}$, and consider the hermitian
Hamiltonian $\hat{H}(\lambda,\delta_{\rm in})$. Suppose that there exists some
particular value $\lambda_{\rm in}\in{\mathbb{R}}$ at which the eigenvalue
spectrum of $\hat{H}(\lambda,\delta_{\rm in})$ contains $M_{\rm in}$ simple
binary777 One may of course analyze also more general situations of multiple
degeneracies, but this is beyond the scope of the present paper.  crossings
($1\leq M_{\rm in}\leq N/2$). Meaning that

$\displaystyle\;\;\;\;\;E_{1}(\lambda_{\rm in},\delta_{\rm in})$
$\displaystyle=$ $\displaystyle\;\;E_{2}(\lambda_{\rm in},\delta_{\rm
in})\hskip 14.22636pt;$ $\displaystyle\;\;\;\;\;E_{3}(\lambda_{\rm
in},\delta_{\rm in})$ $\displaystyle=$ $\displaystyle\;\;E_{4}(\lambda_{\rm
in},\delta_{\rm in})\hskip 14.22636pt;$ $\displaystyle\;\;\bm{\vdots}$
$\displaystyle E_{2M_{\rm in}-1}(\lambda_{\rm in},\delta_{\rm in})$
$\displaystyle=$ $\displaystyle E_{2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})\hskip 7.11317pt;$

where the twice degenerate levels $E_{1}(\lambda_{\rm in},\delta_{\rm
in}),\ldots,E_{2M_{\rm in}-1}(\lambda_{\rm in},\delta_{\rm in})$ are all
distinct, and satisfy also

$\displaystyle\;\;\;\;\;\partial_{\lambda}\,E_{1}(\lambda,\delta_{\rm
in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}$ $\displaystyle\neq$
$\displaystyle\;\;\partial_{\lambda}\,E_{2}(\lambda,\delta_{\rm
in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}\hskip 14.22636pt;$
$\displaystyle\;\;\;\;\;\partial_{\lambda}\,E_{3}(\lambda,\delta_{\rm
in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}$ $\displaystyle\neq$
$\displaystyle\;\;\partial_{\lambda}\,E_{4}(\lambda,\delta_{\rm
in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}\hskip 14.22636pt;$
$\displaystyle\;\;\bm{\vdots}$ $\displaystyle\partial_{\lambda}\,E_{2M_{\rm
in}-1}(\lambda,\delta_{\rm in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}$
$\displaystyle\neq$ $\displaystyle\partial_{\lambda}\,E_{2M_{\rm
in}}(\lambda,\delta_{\rm in})\,\Bigr{|}_{\lambda=\lambda_{\rm in}}\hskip
7.11317pt.$

Conditions (A) say that each of the listed ”simple binary” crossings
corresponds to an intersection of two $\lambda$-dependent eigenvalue lines
with nonequal slopes. All the remaining eigenvalues $E_{j>2M_{\rm
in}}(\lambda_{\rm in},\delta_{\rm in})$ are assumed to be nondegenerate. Figs.
1, 2, 3, 7 in the main text illustrate neatly the presence of the just
discussed simple binary crossings in the case of our toy model Hamiltonian at
$\delta_{\rm in}=0$. In fact, Figs. 1, 2, 3, 7 depict even several distinct
occurrences of $\lambda_{\rm in}$ together with their pertinent values of
$M_{\rm in}$ (one may actually observe that $M_{\rm in}\in\\{1,2,4\\}$ in
these plots).

Let us explore now what happens with a particular multiplet of simple binary
crossings $(\delta_{\rm in},\lambda_{\rm in},M_{\rm in})$ once
$\delta\in{\mathbb{R}}$ is set to depart slightly from $\delta_{\rm in}$, and
once $\lambda$ is set to deviate slightly from $\lambda_{\rm in}$ while being
allowed to penetrate into the complex plane. As a matter of fact, each
crossing $\kappa\in\\{1,2,\bm{\cdots},M_{\rm in}\\}$ survives inside the
complex $\lambda$-plane in the form of a binary EP, which moves with $\delta$
along a certain well defined trajectory $(\delta,\lambda_{\kappa}(\delta))$.
Generally speaking, the resulting trajectories $\lambda_{\kappa}(\delta)$ will
be $\kappa$-dependent. However, eventual symmetries of
$\hat{H}(\lambda,\delta)$ may also cause (some of) these trajectories to be
exactly identical. Under these more peculiar circumstances, our $M_{\rm in}$
binary crossings can be classified into subgroups (clusters), such that
$\lambda_{\kappa}(\delta)$ is the same within each subgroup (cluster).
Different clusters will be hereafter labeled by index $k$.

Consider now a particular $k$-th cluster of $M$ binary EPs $(1\leq M\leq
M_{\rm in})$. As explained in the previous paragraph, this cluster of $M$ EPs
(whose elements we are going to label by index $m\in\\{1,2,\bm{\cdots},M\\}$)
emanates from a subset of simple binary hermitian crossings $(\delta_{\rm
in},\lambda_{\rm in},M_{\rm in})$, and $1\leq M\leq M_{\rm in}$. All the
mentioned $M$ EPs are associated with the same complex $\lambda$-trajectory,
$\lambda_{k}(\delta)$. At $\delta=\delta_{\rm in}$, one has

$\lambda_{k}(\delta_{\rm in})\;=\;\lambda_{\rm in}\hskip 14.22636pt;$ (65)

and

$\displaystyle\;\;\;\;\;E_{1}(\lambda_{\rm in},\delta_{\rm in})$
$\displaystyle=$ $\displaystyle\;\;E_{2}(\lambda_{\rm in},\delta_{\rm
in})\hskip 14.22636pt;$ $\displaystyle\;\;\;\;\;E_{3}(\lambda_{\rm
in},\delta_{\rm in})$ $\displaystyle=$ $\displaystyle\;\;E_{4}(\lambda_{\rm
in},\delta_{\rm in})\hskip 14.22636pt;$ $\displaystyle\;\;\bm{\vdots}$
$\displaystyle E_{2M-1}(\lambda_{\rm in},\delta_{\rm in})$ $\displaystyle=$
$\displaystyle E_{2M}(\lambda_{\rm in},\delta_{\rm in})\hskip 14.22636pt.$

We have conveniently adopted here the same kind of notation as above in (A).

Let the orthonormalized eigenvectors corresponding to $E_{1}(\lambda_{\rm
in},\delta_{\rm in})$, $E_{2}(\lambda_{\rm in},\delta_{\rm in})$,
$E_{3}(\lambda_{\rm in},\delta_{\rm in})$ etc. be denoted by symbols
$|v_{1}(\lambda_{\rm in},\delta_{\rm in})\rangle$, $|v_{2}(\lambda_{\rm
in},\delta_{\rm in})\rangle$, $|v_{3}(\lambda_{\rm in},\delta_{\rm
in})\rangle$, etc. Note that we use here the standard ket-notation, since
$\hat{H}(\lambda_{\rm in},\delta_{\rm in})$ is hermitian (real symmetric) and
thus the conventional definition of the scalar product applies. Since
$E_{1}(\lambda_{\rm in},\delta_{\rm in})=E_{2}(\lambda_{\rm in},\delta_{\rm
in})$, the sought initial condition for the $m=1$ EP must inevitably look as
follows:

$\tilde{E}_{1}^{\delta_{\rm in}}\;=\;E_{1}(\lambda_{\rm in},\delta_{\rm
in})\hskip 14.22636pt;$ (67)

and

$\displaystyle|\tilde{c}_{1}^{\delta_{\rm in}})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,\Bigl{(}|v_{1}(\lambda_{\rm in},\delta_{\rm
in})\rangle\,+\,\sigma_{1}\,i\,|v_{2}(\lambda_{\rm in},\delta_{\rm
in})\rangle\Bigr{)}\hskip 14.22636pt;$ (68)
$\displaystyle|\tilde{b}_{1}^{\delta_{\rm in}})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,\Bigl{(}|v_{1}(\lambda_{\rm in},\delta_{\rm
in})\rangle\,-\,\sigma_{1}\,i\,|v_{2}(\lambda_{\rm in},\delta_{\rm
in})\rangle\Bigr{)}\hskip 14.22636pt.$ (69)

In (68)-(69), the sign factor $\sigma_{1}\in\\{-1,+1\\}$. Similarly for all
the other EPs $m=2,3,\bm{\cdots},M$. Written down explicitly, we set

$\tilde{E}_{m}^{\delta_{\rm in}}\;=\;E_{2m-1}(\lambda_{\rm in},\delta_{\rm
in})\hskip 14.22636pt;$ (70)

and

$\displaystyle|\tilde{c}_{m}^{\delta_{\rm in}})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,\Bigl{(}|v_{2m-1}(\lambda_{\rm
in},\delta_{\rm in})\rangle\,+\,\sigma_{m}\,i\,|v_{2m}(\lambda_{\rm
in},\delta_{\rm in})\rangle\Bigr{)}\hskip 14.22636pt;$ (71)
$\displaystyle|\tilde{b}_{m}^{\delta_{\rm in}})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,\Bigl{(}|v_{2m-1}(\lambda_{\rm
in},\delta_{\rm in})\rangle\,-\,\sigma_{m}\,i\,|v_{2m}(\lambda_{\rm
in},\delta_{\rm in})\rangle\Bigr{)}\hskip 14.22636pt;$ (72)

where $1\leq m\leq M$ and $\sigma_{m}\in\\{-1,+1\\}$. An assignment of the
sign factors $(\sigma_{1},\sigma_{2},\bm{\cdots},\sigma_{M})$ in (71)-(72)
must be performed in such a consistent way that the velocity

$\bm{\dot{}}{\lambda}(\delta_{\rm
in})\;=\;-\,\frac{(\tilde{c}_{m}^{\delta_{\rm
in}}|\partial_{\delta}\,\hat{H}(\lambda_{\rm in},\delta_{\rm
in})|\tilde{c}_{m}^{\delta_{\rm in}})}{(\tilde{c}_{m}^{\delta_{\rm
in}}|\partial_{\lambda}\,\hat{H}(\lambda_{\rm in},\delta_{\rm
in})|\tilde{c}_{m}^{\delta_{\rm in}})}$ (73)

predicted by equation (24) comes out as being independent of $m$. We shall
return to the sign factors $(\sigma_{1},\sigma_{2},\bm{\cdots},\sigma_{M})$
below (see the item (ii) in the last paragraph).

We need to specify also the ICS for all the ordinary non-EP eigenstates of
$\hat{H}(\lambda_{k}(\delta_{\rm in}),\delta_{\rm in})$. This task is
straightforward in the case of non-degenerate energy levels $E_{j+2M_{\rm
in}}(\lambda_{\rm in},\delta_{\rm in})$ (where $1\leq j\leq N-2M_{\rm in}$).
One sets obviously

$E_{j}^{\delta_{\rm in}}\;=\;E_{j+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})\hskip 14.22636pt;\hskip 14.22636pt1\leq j\leq N-2M_{\rm in}$ (74)

and

$|c_{j}^{\delta_{\rm in}})\;=\;|v_{j+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})\rangle\hskip 14.22636pt;\hskip 14.22636pt1\leq j\leq N-2M_{\rm in}$ (75)

where $|v_{j+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm in})\rangle$ stands of
course for the unit normalized eigenvector of $\hat{H}(\lambda_{\rm
in},\delta_{\rm in})$ associated with level $E_{j+2M_{\rm in}}(\lambda_{\rm
in},\delta_{\rm in})$.

The situation becomes somewhat more delicate in the case of the doubly
degenerate energy eigenvalues listed in (A) but not included in the cluster
(A), namely, in the case of levels

$\displaystyle\;\;E_{2M+1}(\lambda_{\rm in},\delta_{\rm in})$ $\displaystyle=$
$\displaystyle E_{2M+2}(\lambda_{\rm in},\delta_{\rm in})\hskip 14.22636pt;$
$\displaystyle\;\;E_{2M+3}(\lambda_{\rm in},\delta_{\rm in})$ $\displaystyle=$
$\displaystyle E_{2M+4}(\lambda_{\rm in},\delta_{\rm in})\hskip 14.22636pt;$
$\displaystyle\;\;\bm{\vdots}$ $\displaystyle E_{2M_{\rm in}-1}(\lambda_{\rm
in},\delta_{\rm in})$ $\displaystyle=$ $\displaystyle E_{2M_{\rm
in}}(\lambda_{\rm in},\delta_{\rm in})\hskip 14.22636pt\;\,.$

Consider any given doubly degenerate eigenvalue

$E_{2M+j-N+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})\;=\;E_{2M+j-N+2M_{\rm in}+1}(\lambda_{\rm in},\delta_{\rm in})$ (76)

where $N-2M_{\rm in}+1\leq j\leq N-2M-1$. Let the two pertinent unit
normalized orthonormal eigenvectors be

$\displaystyle|v^{(1)}\rangle$ $\displaystyle\equiv$
$\displaystyle|v_{2M+j-N+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})\rangle\hskip 14.22636pt\;\;\;\,;$ (77) $\displaystyle|v^{(2)}\rangle$
$\displaystyle\equiv$ $\displaystyle|v_{2M+j-N+2M_{\rm in}+1}(\lambda_{\rm
in},\delta_{\rm in})\rangle\hskip 14.22636pt.$ (78)

We set of course

$E_{2M+j-N+2M_{\rm in}}^{\delta_{\rm in}}\;=\;E_{2M+j-N+2M_{\rm
in}}(\lambda_{\rm in},\delta_{\rm in})\;=\;E_{2M+j-N+2M_{\rm
in}+1}^{\delta_{\rm in}}\hskip 14.22636pt;$ (79)

much as in (74). Yet an assignment of the corresponding non-EP eigenvectors

$|c_{2M+j-N+2M_{\rm in}}^{\delta_{\rm in}})\hskip 14.22636pt{\rm and}\hskip
14.22636pt|c_{2M+j-N+2M_{\rm in}+1}^{\delta_{\rm in}})$ (80)

needs a bit more care. Clearly, entities (80) must be built up as
$c$-orthonormalized linear combinations of the two eigenstates (77)-(78). In
addition, however, one must ensure that the two sought non-EP eigenvectors
(80) are not mutually coupled by the Hamiltonian $\delta$-derivative (7),
i.e., by the operator

$\hat{V}(\delta_{\rm in})\;=\;\partial_{\lambda}\,\hat{H}(\lambda_{\rm
in},\delta_{\rm in})\,\bm{\dot{}}{\lambda}(\delta_{\rm
in})\;+\;\partial_{\delta}\,\hat{H}(\lambda_{\rm in},\delta_{\rm in})\hskip
14.22636pt.$ (81)

Indeed, the just imposed extra requirement of

$(c_{2M+j-N+2M_{\rm in}}^{\delta_{\rm in}}|\hat{V}(\delta_{\rm
in})|c_{2M+j-N+2M_{\rm in}+1}^{\delta_{\rm in}})\;=\;0$ (82)

is indispensable, since it guarantees that our EOM (45) does not possess a
singularity at $\delta=\delta_{\rm in}$. Hence an appropriate kind of
regularization or rectification must be implemented here. In fact, an explicit
construction of the two non-EP eigenvectors (80) is conceptually
straightforward. Namely, we diagonalize888 The matrix (83) is surely
diagonalizable. Since if it was non-diagonalizable, then the just investigated
crossing of eigenvalues $E_{2M+j-N+2M_{\rm in}}(\lambda_{\rm in},\delta_{\rm
in})=E_{2M+j-N+2M_{\rm in}+1}(\lambda_{\rm in},\delta_{\rm in})$ would bring
an additional $(M+1)$-th EP into the list (A), contrary to our starting
assumption.  the 2-by-2 matrix

$\left(\matrix{\langle v^{(1)}|\hat{V}(\delta_{\rm in})|v^{(1)}\rangle&\langle
v^{(1)}|\hat{V}(\delta_{\rm in})|v^{(2)}\rangle\cr\langle
v^{(2)}|\hat{V}(\delta_{\rm in})|v^{(1)}\rangle&\langle
v^{(2)}|\hat{V}(\delta_{\rm in})|v^{(2)}\rangle}\right)\hskip 14.22636pt;$
(83)

and access in this way the two associated eigenvectors

$\vec{w}^{(1)}\;=\;\left(\matrix{w_{1}^{(1)}\cr w_{2}^{(1)}}\right)\hskip
14.22636pt,\hskip 14.22636pt\vec{w}^{(2)}\;=\;\left(\matrix{w_{1}^{(2)}\cr
w_{2}^{(2)}}\right)\hskip 14.22636pt.$ (84)

Subsequently we set

$\displaystyle|c_{2M+j-N+2M_{\rm in}\phantom{+1}}^{\delta_{\rm in}})$
$\displaystyle=$ $\displaystyle
w_{1}^{(1)}\,|v^{(1)}\rangle\;+\;w_{2}^{(1)}\,|v^{(2)}\rangle\hskip
14.22636pt;$ (85) $\displaystyle|c_{2M+j-N+2M_{\rm in}+1}^{\delta_{\rm in}})$
$\displaystyle=$ $\displaystyle
w_{1}^{(2)}\,|v^{(1)}\rangle\;+\;w_{2}^{(2)}\,|v^{(2)}\rangle\hskip
14.22636pt;$ (86)

while tacitly implementing the $c$-normalization.

What remains to be done is to supply the values of $f_{m}^{\delta_{\rm in}}$.
Equation (10) implies

$f_{m}^{\delta_{\rm in}}\;=\;0\hskip 14.22636pt;\hskip 14.22636pt1\leq m\leq
M$ (87)

valid simply because $|\tilde{b}_{m}^{\delta_{\rm in}})$ is an eigenvector of
$\hat{H}(\lambda_{k}(\delta_{\rm in}),\delta_{\rm in})$ with an eigenvalue
$\tilde{E}_{m}^{\delta_{\rm in}}$.

Summarizing, in this Appendix A we have described in a self contained fashion
how to specify adequately the ICS for the seven fundamental entities (19). The
resulting ICS are given above in equations (65), (70), (71), (72), (74), (75),
(79), (85), (86), and (87). The just mentioned ICS must satisfy by
construction the basic eigenvalue and eigenvector properties (8), (9), (10),
(11)-(16), (17) listed in Section 2 of the main text, this may serve as an
useful consistency check.

Let us finally mention two important questions which still need to be
addressed:

  * (i)

A nontrivial puzzle arises on how a given multiplet of simple binary crossings
$(\delta_{\rm in},\lambda_{\rm in},M_{\rm in}\geq 2)$ should be split into
specific subgroups (clusters) characterized by the same $\lambda_{k}(\delta)$.

  * (ii)

Another nontrivial puzzle concerns consistent choice of the $M$ sign factors
$(\sigma_{1},\sigma_{2},\bm{\cdots},\sigma_{M})$ for a given $k$-th cluster in
equations (71)-(72), see the above discussion of requirement (73).

Both puzzles (i)-(ii) are resolved correctly iff an explicit solution of our
EOM (24), (26), (31), (43), (44), (45), (46), starting from the just discussed
ICS, provides unique outcomes (19) which do possess the basic properties (8),
(9), (10), (11)-(16), (17) of Section 2 for all values of $\delta$ considered
in the propagation. On the other hand, any inconsistency detected during the
propagation of our EOM, manifested e.g. by violation of any from the
properties (8), (9), (10), (11)-(16), (17), would inevitably imply an
incorrect resolution of one or both of the aforementioned issues (i)-(ii).
Hence the two puzzles (i)-(ii) can be uniquely resolved simply on the trial-
and-error basis, even in such situations when direct answer to (i)-(ii) is not
a priori obvious.