Perturbation Theory Near Degenerate Exceptional Points

Abstract

In an overall framework of quantum mechanics of unitary systems a rather sophisticated new version of perturbation theory is developed and described. The motivation of such an extension of the list of the currently available perturbation-approximation recipes was four-fold: (1) its need results from the quick growth of interest in quantum systems exhibiting parity-time symmetry ($\mathcal{PT}$-symmetry) and its generalizations; (2) in the context of physics, the necessity of a thorough update of perturbation theory became clear immediately after the identification of a class of quantum phase transitions with the non-Hermitian spectral degeneracies at the Kato’s exceptional points (EP); (3) in the dedicated literature, the EPs are only being studied in the special scenarios characterized by the spectral geometric multiplicity L equal to one; (4) apparently, one of the decisive reasons may be seen in the complicated nature of mathematics behind the $L\ge 2$ constructions. In our present paper we show how to overcome the latter, purely technical obstacle. The temporarily forgotten class of the $L>1$ models is shown accessible to a feasible perturbation-approximation analysis. In particular, an emergence of a counterintuitive connection between the value of L, the structure of the matrix elements of perturbations, and the possible loss of the stability and unitarity of the processes of the unfolding of the singularities is given a detailed explanation.

1. Introduction

The Bender’s and Boettcher’s [1] idea of replacement of Hermiticity $H=H^{†}$ by parity-time symmetry ($\mathcal{PT}$-symmetry) $H\mathcal{PT}=\mathcal{PT}H$ of a Hamiltonian responsible for unitary evolution opened, after an appropriate mathematical completion [2,3,4,5] of the theory, a way towards the building of quantum models exhibiting non-Hermitian degeneracies [6] alias exceptional points (EPs, [7]). For a fairly realistic illustrative example of possible applications opening multiple new horizons in phenomenology we could recall, e.g., the well-known phenomenon of Bose–Einstein condensation is a schematic simplification described by the non-Hermitian but $\mathcal{PT}$-symmetric three-parametric Hamiltonian

$$H^{\left(GGKN\right)}(\gamma ,v,c)=-\mathrm{i}\gamma \left(a_1^{†}{a}_1-a_2^{†}{a}_2\right)+v\left(a_1^{†}{a}_2+a_2^{†}{a}_1\right)+\frac{c}{2}{\left(a_1^{†}{a}_1-a_2^{†}{a}_2\right)}^2.$$

(1)

This Hamiltonian represents an interesting analytic-continuation modification of the conventional Hermitian Bose–Hubbard Hamiltonian [8,9,10]. In this form the model was recently paid detailed attention in Ref. [11]. A consequent application of multiple, often fairly sophisticated forms of perturbation theory has been shown there to lead to surprising results. In particular, the behavior of the bound and resonant states of the system was found to lead to the new and unexpected phenomena in the dynamical regime characterized by the small coupling constant c. The authors of Ref. [11] emphasized that new physics may be expected to emerge precisely in the vicinity of the EP-related dynamical singularities.

These phenomena (simulating not necessarily just the Bose–Einstein condensation of course) were analyzed, in [11], using several ad hoc, not entirely standard perturbation techniques. The role of an unperturbed Hamiltonian $H_0$ was assigned, typically, to the extreme EP limits of H. Unfortunately, only too often the perturbed energies appeared to be complex as a consequence. In other words, the systems exhibiting $\mathcal{PT}$-symmetry seemed to favor the spontaneous breakdown of this symmetry near EPs.

In the light of similar results one immediately must ask the question whether such a “wild behavior” of the EP-related quantum systems is generic. Indeed, an affirmative answer is often encountered in the studies by mathematicians (see, e.g., [12]). A hidden reason is that they usually tacitly keep in mind just the “effective theory” and/or the so-called “open quantum system” dynamical scenario [13].

In the more restrictive context of the unitary quantum mechanics the situation is different: the “wild behavior” of systems is usually not generic there (cf. also the recent explanatory commentary on the sources of possible misunderstandings in [14]). Several non-numerical illustrative models may be found in [15] where typically, the $\mathcal{PT}$-symmetric Bose–Hubbard model of Equation (1) (which behaves as unstable near its EP singularities [11])) has been replaced by its “softly perturbed” alternative in which in an arbitrarily small vicinity of its EP singularity, the system remains stable and unitary under admissible perturbations.

The physics of stability covered by paper [15] can be perceived as one of the main sources of inspiration of our present study. We intend to replace here the very specific model of Equation (1) (in which the geometric multiplicity L of all of its EP-related degeneracies was always equal to one) by a broader class of quantum systems. In a way motivated by the idea of a highly desirable extension of the currently available menu of the tractable and eligible dynamical scenarios beyond their $L=1$ subclass, we will turn attention here to the EP-related degeneracies of the larger, nontrivial geometric multiplicities $L\ge 2$. We will reveal that such a study opens new horizons not only in phenomenology (where the influence of perturbations becomes strongly dependent on the detailed structure of the non-Hermitian degeneracy) but also in mathematics (where a rich menu of physical consequences will be shown reflected by an unexpected adaptability of the geometry of the Hilbert space to the detailed structure of the perturbation).

The presentation of our results will be organized as follows. First, in Section 2 we will recall a typical quantum system (viz., a version of the non-Hermitian Bose–Hubbard multi-bosonic model) in which the EP degeneracies play a decisive phenomenological role. We will explain that although the model itself only exhibits the maximal-order $L=1$ EP degeneracies, such an option represents, from the purely formal point of view, just one of the eligible dynamical scenarios. In Appendix A a full classification of the EPs is presented therefore, showing, i.a., that the number of the “anomalous” EPs of our present interest with $L\ge 2$ exhibits an almost exponential growth at the larger matrix dimension N.

The goals of our considerations are subsequently explained in Section 3. For the sake of brevity, we just pick up the first nontrivial case with $L=2$, and we emphasize that even in such a case the basic features of an appropriate adaptation of perturbation theory may be explained, exhibiting also, not quite expectedly, a survival of the fairly user-friendly mathematical structure.

In order to make our message self-contained, the known form of the EP-related perturbation formalism restricted to $L=1$ is reviewed in Appendix B. On this background, in a way based on a not too dissimilar constructive strategy, our present main $L\ge 2$ results are then presented and described in Section 4. We emphasize there the existence of the phenomenological as well as mathematical subtleties of the large-L models. We show that in our generalized, degenerate-perturbation-theory formalism a key role is played by an interplay between its formal mathematical background (viz, the non-Hermiticity of the Hamiltonians) and its phenomenological aspects (typically, the knowledge of L must be complemented by an explicit knowledge of the partitioning of Schrödinger equation).

In Section 5, all these aspects of the $L>1$ perturbation theory are summarized and illustrated by a detailed description of the characteristic, not always expected features of the leading-order approximations. Several related applicability aspects of our present degenerate-perturbation-theory formalism are finally discussed in Section 6 and in two Appendices. We point out there that some of the features of the $L>1$ theory (e.g., a qualitative, fairly counterintuitive clarification of the concept of the smallness of perturbations) may be treated as not too different from their $L=1$ predecessors. At the same time, a wealth of new formal challenges is emphasized to emerge, in particular, in the analyses of the role of perturbations in the specific quantum systems which are required unitary.

2. Exceptional Points

2.1. Bose–Hubbard Model and Exceptional Points of Geometric Multiplicity One, $L=1$

Illustrative $\mathcal{PT}$-symmetric Bose–Hubbard Hamiltonian operator (1) commutes with the number operator

$$\widehat{N}=a_1^{†}{a}_1+a_2^{†}{a}_2.$$

(2)

This means that the number of bosons N is conserved so that after its choice the Hamiltonian may be represented by a finite-dimensional non-Hermitian K by K matrix $H^{\left(GGKN\right)}(K,\gamma ,v,c)$ with $K=N+1$ (see its explicit construction in [11]). Once we fix the units (such that $v=1$) and once we set $c=0$ (preserving, for the sake of simplicity, just the first two components of the Hamiltonian), the resulting one-parametric family of Hamiltonian matrices $H^{\left(GGKN\right)}(K,\gamma )$ can be assigned the closed-form energy spectra

$$E_n^{\left(GGKN\right)}(K,\gamma )={(1-{\gamma }^2)}^{1/2}(1-K+2n),n=0,1,\dots ,K-1.$$

(3)

These energies remain real and non-degenerate (i.e., observable) if and only if ${\gamma }^2<1$.

In loc. cit. it has also been proved that the two interval-boundary values of $\gamma =\pm 1$ are, in the terminology of the Kato’s mathematical monograph [7], exceptional points (EPs, ${\gamma }_{+}^{\left(EP\right)}=1$ and ${\gamma }_{-}^{\left(EP\right)}=-1$). More precisely, one should speak about the very special EPs of maximal order (i.e., of order K, abbreviated as EPK). The latter observation may be given a more general, model-independent linear-algebraic background via relation

$$H^{\left(K\right)}({\gamma }^{\left(EPK\right)})Q^{\left(EPK\right)}=Q^{\left(EPK\right)}{J}^{\left(K\right)}\left(\eta \right)$$

(4)

where

$$\eta =\underset{\gamma \to {\gamma }^{\left(EPK\right)}}{lim}{E}_n(K,\gamma ),n=0,1,\dots ,K-1$$

(5)

is the limiting degenerate energy. Relation (4) contains the so-called transition matrix $Q^{\left(EPK\right)}$ and the non-diagonal, EPK-related canonical-representation Jordan-block matrix

$$J^{\left(K\right)}\left(\eta \right)=\left(\begin{array}{ccccc}\eta & 1& 0& \dots & 0\\ 0& \eta & 1& \ddots & ⋮\\ 0& 0& \eta & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & 1\\ 0& \dots & 0& 0& \eta \end{array}\right).$$

(6)

As a certain limiting analogue of the conventional set of eigenvectors the transition matrix is obtainable via the solution of the EPK-related analogue (4) of conventional Schrödinger equation. For our illustrative example $H^{\left(GGKN\right)}(K,\gamma )$, in particular, all of the K- and ${\gamma }^{\left(EPK\right)}$-dependent explicit, closed forms of solutions $Q^{\left(EPK\right)}$ remain non-numerical and may be found constructed in dedicated paper [15].

2.2. Generic Non-Hermitian Degeneracies with Geometric Multiplicities Larger than One, $L>1$

In the common model-building scenarios the $\gamma$-dependence of Hamiltonians $H^{\left(K\right)}\left(\gamma \right)$ is analytic. Under this assumption the exceptional points of maximal order as discussed in preceding subsection represent just one of several possible realizations of a non-Hermitian degeneracy (NHD) with its characteristic EP-related confluence of eigenvalues (5). Besides the maximal, EPK-related complete confluence of eigenvectors as described in preceding subsection we may encounter, in general, multiple other, incomplete confluences of eigenvectors

$$\underset{\gamma \to {\gamma }^{\left(NHD\right)}}{lim}|{\psi }_{m_{k\left[j\right]}}\left(\gamma \right)\rangle =|{\chi }_j\rangle ,k\left[j\right]=1,2,\dots ,M_j,M_j\ge 2,j=1,2,\dots ,L.$$

(7)

Here we have $M_1+M_2+\dots +M_L=K$ where L is called the geometric multiplicity of the EP degeneracy [7] (see also Appendix A for more details).

Every EP instant ${\gamma }^{\left(NHD\right)}$ may be characterized not only by the overall Hilbert-space dimension K and by the number $L\ge 1$ of linearly independent $\eta$-related states $|{\chi }_j\rangle$ of Equation (7) but also by a suitably ordered L-plet of the related subspace dimensions $M_j$. Thus, in the present $L\ge 2$ extension of preceding subsection the fully non-diagonal Jordan block of Equation (6) must be replaced by the more general block-diagonal canonical representation of the Hamiltonian,

$${\mathcal{J}}^{\left(K\right)}\left(\eta \right)=\left(\begin{array}{cccc}{J}^{\left(M_1\right)}\left(\eta \right)& 0& \dots & 0\\ 0& J^{\left(M_2\right)}\left(\eta \right)& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \dots & 0& J^{\left(M_L\right)}\left(\eta \right)\end{array}\right)=\underset{j=1}{\overset{L}{⨁}}{J}^{\left(M_j\right)}\left(\eta \right).$$

(8)

In parallel, EPK relation (4) must be replaced by its generalization

$$H^{\left(K\right)}({\gamma }^{\left(NHD\right)})Q^{\left(NHD\right)}=Q^{\left(NHD\right)}{\mathcal{J}}^{\left(K\right)}\left(\eta \right).$$

(9)

Naturally, the direct-sum structure of ${\mathcal{J}}^{\left(K\right)}\left(\eta \right)$ becomes reflected by a partitioned-matrix structure of transition matrices $Q^{\left(NHD\right)}$ which are, in general, not block-diagonal of course.

3. Unitary Processes of Collapse at L = 2

A priori one may expect that the existence of anisotropy of the Hilbert space as realized, in Equation (A18) at $L=1$, by an elementary rescaling $B\left(\lambda \right)$ of the basis will also exist at any larger geometric multiplicity $L>1$, i.e., in the unitary quantum systems with the more complicated structure of the NHD limiting alias quantum phase transition.

3.1. Quantum Physics behind “Degenerate Degeneracies” with $L=2$

Hypothetically, the unitary evolution of any $\mathcal{PT}$-symmetric quantum system moving towards a hiddenly Hermitian EP degeneracy with geometric multiplicity two can be perceived as generated by a suitable diagonalizable Hamiltonian $H^{\left(K\right)}\left(\gamma \right)$ with real spectrum [3]. What is only necessary is that its (perhaps, properly renumbered) eigenvectors $|{\mathsf{\Phi }}_j\left(\gamma \right)\rangle$ obey the $L=2$ EP-degeneracy rule

$$\underset{\gamma \to {\gamma }^{\left(EP\right)}}{lim}|{\mathsf{\Phi }}_m\left(\gamma \right)\rangle =|{\chi }_a\rangle ,m=0,1,\dots ,M-1,$$

(10)

$$\underset{\gamma \to {\gamma }^{\left(EP\right)}}{lim}|{\mathsf{\Phi }}_{M+n}\left(\gamma \right)\rangle =|{\chi }_b\rangle ,n=0,1,\dots ,N-1.$$

(11)

The two limiting eigenvectors $|{\chi }_a\rangle$ and $|{\chi }_b\rangle$ are, by our assumption, linearly independent so that partitioned matrix (8), i.e., at $L=2$ and $\eta =0$, matrix

$${\mathcal{J}}^{(M⨁N)}\left(0\right)=\left[\begin{array}{cc}{J}^{\left(M\right)}\left(0\right)& 0\\ 0& J^{\left(N\right)}\left(0\right)\end{array}\right]$$

(12)

will represent the canonical form of our Hamiltonian in the EP limit. The corresponding transition matrix $Q^{(M⨁N)}$ can be then obtained by the solution of the limiting version

$$H^{\left(K\right)}({\gamma }^{\left(EP\right)})Q^{(M⨁N)}=Q^{(M⨁N)}{\mathcal{J}}^{(M⨁N)}\left(0\right)$$

(13)

of the initial Schrödinger equation. A few exactly solvable samples of the latter $L=2$ degeneracy process $\gamma \to {\gamma }^{\left(EP\right)}$ can be found, e.g., in our recent paper [16].

In all the similar dynamical scenarios one can always find parallels with their simpler $L=1$ predecessors. In particular, a return to the situation before the collapse as sampled, at $L=1$, by Equation (A4) below, can be also given the following analogous form

$${\left[Q^{\left(\right(M⨁N)}\right]}^{-1}{H}^{\left(K\right)}\left(\gamma \right)Q^{(M⨁N)}=J^{(M⨁N)}\left(0\right)+\lambda V^{\left(K\right)}\left(\gamma \right),\gamma \ne {\gamma }^{\left(EP\right)}$$

(14)

of a “unitarity-compatible” perturbation-theoretic reinterpretation in which the perturbation $\lambda V^{\left(K\right)}\left(\gamma \right)$ is fully determined by the input matrix Hamiltonian $H^{\left(K\right)}\left(\gamma \right)$.

3.2. Unfoldings of Degeneracies under Random Perturbations at $L=2$

In a close parallel to the $L=1$ Schrödinger’s bound-state problem (A6) let us now start the study of its $L>1$ generalizations by considering the first nontrivial choice of degeneracy with the geometric multiplicity $L=2$. In our present notation the corresponding Schrödinger equation then reads

$$\left[{\mathcal{J}}^{(M⨁N)}\left(0\right)+\lambda V^{\left(K\right)}\right]|\mathsf{\Psi }\rangle =\epsilon |\mathsf{\Psi }\rangle .$$

(15)

Without any loss of generality, we set again $\eta =0$. After such a choice all the eigenvalues $\epsilon =\epsilon \left(\lambda \right)$ will remain small, and they will vanish in the formal unperturbed-system limit $\lambda \to 0$.

In a way paralleling Equation (12), any given matrix of perturbations has to be partitioned as well,

$$V_{}^{\left(K\right)}=\left[\begin{array}{cc}{V}_{}^{(M,M)}& V_{}^{(M,N)}\\ V_{}^{(N,M)}& V_{}^{(N,N)}\end{array}\right].$$

(16)

Since the four submatrices have dimensions indicated by the superscripts, we shall assume, at the beginning at least, that all the individual matrix elements of the perturbation matrix $V^{\left(K\right)}$ remain bounded at small $\lambda$. This means that the Hamiltonian is dominated by its unperturbed part (12) so that also the whole perturbed $L=2$ Schrödinger Equation (15) has to be partitioned. In order to simplify the notation we shall write

$$|\mathsf{\Psi }\rangle =\left(\begin{array}{c}|{\psi }^{\left(a\right)}\succ \\ |{\psi }^{\left(b\right)}\succ \end{array}\right),|{\psi }^{\left(a\right)}\succ =\left(\begin{array}{c}{\psi }_1^{\left(a\right)}\\ {\psi }_2^{\left(a\right)}\\ ⋮\\ {\psi }_M^{\left(a\right)}\end{array}\right),|{\psi }^{\left(b\right)}\succ =\left(\begin{array}{c}{\psi }_1^{\left(b\right)}\\ {\psi }_2^{\left(b\right)}\\ ⋮\\ {\psi }_N^{\left(b\right)}\end{array}\right)$$

(17)

using the curly-ket symbols for subvectors. This will enable us to proceed in a partial parallel with the widely studied non-degenerate-EP cases where one has $L=1$ (see Refs. [17,18] or Appendix B below for a compact review).

4. Perturbation Theory at L = 2

4.1. The Recent Change of the Unitary-Evolution Paradigm

From the historical point of view the use of EPs in physics has not been immediate. Only during the last circa 20 years one notices a perceivable increase of the relevance of the concept in various branches of theoretical as well as experimental physics. Various innovative EP applications emerged ranging from the analyses of resonances in classical mechanics [19] and of the so called non-Hermitian degeneracies in classical optics [6] up to the studies of the wealth of phenomena in quantum physics of open quantum systems [13,20] or, finally, even of the closed, stable and unitarily evolving quantum systems [1,2].

All these developments contributed to the motivation of our present study. For the sake of definiteness, we restricted our attention to the framework of quantum physics, unitary or non-unitary. In this setting the traditional role of the EPs ${\gamma }^{\left(EP\right)}$ has always been two-fold. First, in the context of mathematics, the conventional analyticity assumptions about Hamiltonians

$$H\left(\gamma \right)=H\left({\gamma }_0\right)+(\gamma -{\gamma }_0)H^{\left(1\right)}+{(\gamma -{\gamma }_0)}^2{H}^{\left(2\right)}+\dots ,$$

(18)

and the conventional power-series ansatz for energies

$$E_n\left(\gamma \right)=E_n\left({\gamma }_0\right)+(\gamma -{\gamma }_0)E_n^{\left(1\right)}+{(\gamma -{\gamma }_0)}^2{E}_n^{\left(2\right)}+\dots$$

(19)

(etc.) gave birth to the so-called Rayleigh–Schrödinger perturbation-expansion constructions of the Schrödinger-equation solutions. It has been revealed that the radius of convergence R of these perturbation-series solutions is determined by the position of the nearest EP in the complex plane of the parameter, $R=min|{\gamma }_0-{\gamma }^{\left(EP\right)}|$ [7].

In the other, direct applications of EPs, the localization of singularities ${\gamma }^{\left(EP\right)}$ only played an important traditional role in non-unitary, open quantum systems [13]. An explanation is easy: for any self-adjoint Hamiltonian $H\left(\gamma \right)$ characterizing a closed quantum system, the necessary reality of the parameter ($\mathrm{Im}\gamma =0$) cannot be made compatible with the fact that all of the values of the eligible (i.e., not accumulation-point) EP parameters are complex, $\mathrm{Im}{\gamma }^{\left(EP\right)}\ne 0$.

The traditional paradigm has only been changed recently, after Bender with Boettcher [1] managed to turn attention of physicists’ community to the existence of a broad class of Hamiltonians $H\left(\gamma \right)$ which happen to be non-Hermitian but parity-time symmetric ($\mathcal{PT}$-symmetric) in $\mathcal{K}$. One of the characteristic mathematical features of these Hamiltonians is that despite their non-Hermiticity, their whole spectrum $\left\{E_n\right\}$ may remain strictly real in a suitable real interval $\mathcal{D}$ of the unitarity-compatible parameters $\gamma$ (see, e.g., monograph [4] for more details).

4.2. Rearrangement of Schrödinger Equation

With the two subscripts $j_{\left(a\right)}$ and $j_{\left(b\right)}$ running, in the curly-ket subvectors in (17), from 1 to $N_{\left(a\right)}=M$ and $N_{\left(b\right)}=N$, respectively, we will now only partially fix the norm by setting ${\psi }_1^{\left(a\right)}={\omega }_{\left(a\right)}$ and ${\psi }_1^{\left(b\right)}={\omega }_{\left(b\right)}$ or, in a self-explanatory shorthand, ${\psi }_1^{(a,b)}={\omega }_{(a,b)}$. Next, a parallel to the $L=1$ redefinition (A8) of wave functions will be found in its $L=2$ extension

$$|{\overrightarrow{y}}^{(a,b)}\succ =\left(\begin{array}{c}{y}_1^{(a,b)}\\ y_2^{(a,b)}\\ ⋮\\ y_{N_{(a,b)}-1}^{(a,b)}\\ y_{N_{(a,b)}}^{(a,b)}\end{array}\right)=\left(\begin{array}{c}{\psi }_2^{(a,b)}\\ {\psi }_3^{(a,b)}\\ ⋮\\ {\psi }_{N_{(a,b)}}^{(a,b)}\\ {\mathsf{\Omega }}_{(a,b)}\end{array}\right)$$

(20)

where the two new, temporarily variable elements ${\mathsf{\Omega }}_{\left(a\right)}$ and ${\mathsf{\Omega }}_{\left(b\right)}$ will have to be determined later. In terms of the four auxiliary symbols

$$|e^{(a,b)}\succ =\left(\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right),{\mathsf{\Pi }}^{(a,b)}=\left[\begin{array}{ccccc}0& 0& 0& \dots & 0\\ 1& 0& 0& \ddots & ⋮\\ 0& \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & 1& 0& 0\\ 0& \dots & 0& 1& 0\end{array}\right]$$

(21)

of dimensions $N_{\left(a\right)}=M$ and $N_{\left(b\right)}=N$ we will further decompose

$$|{\psi }^{(a,b)}\succ =|e^{(a,b)}\succ {\omega }_{(a,b)}+{\mathsf{\Pi }}^{(a,b)}|{\overrightarrow{y}}^{(a,b)}\succ .$$

In the next step we introduce the $L=2$ analogue of the unpartitioned $L=1$ vector (A7),

$$|r\rangle =\left(\begin{array}{c}|r^{\left(a\right)}\succ \\ |r^{\left(b\right)}\succ \end{array}\right),|r^{(a,b)}\succ =\left(\begin{array}{c}{r}_1^{(a,b)}\\ r_2^{(a,b)}\\ ⋮\\ r_{N_{(a,b)}}^{(a,b)}\end{array}\right)$$

(22)

with components

$$r_i^{(a,b)}=\epsilon {\omega }_{(a,b)}{\delta }_{i,1}-\lambda V_{i,1}^{(N_{(a,b)},M)}{\omega }_{\left(a\right)}-\lambda V_{i,1}^{(N_{(a,b)},N)}{\omega }_{\left(b\right)}=r_i^{(a,b)}(\epsilon ,\overrightarrow{\omega }).$$

(23)

Treating, temporarily, the two not yet specified quantities ${\mathsf{\Omega }}_{\left(a\right)}$ and ${\mathsf{\Omega }}_{\left(b\right)}$ as adjustable matrix-regularization parameters, and replacing the $L=1$ auxiliary matrix (A9) by its partitioned $L=2$ counterpart

$$\mathcal{A}\left(\epsilon \right)=\left[\begin{array}{cc}A(M,\epsilon )& 0\\ 0& A(N,\epsilon )\end{array}\right]$$

we are just left with the problem of finding a suitable $L=2$ analogue of relation (A11).

A key to the resolution of the puzzle is found in Equation (21) and in its partitioned direct-sum extension

$${\mathsf{\Pi }}^{(M⨁N)}=\left[\begin{array}{cc}{\mathsf{\Pi }}^{\left(a\right)}& 0\\ 0& {\mathsf{\Pi }}^{\left(b\right)}\end{array}\right].$$

Using this symbol, we can now rewrite our homogeneous Schrödinger Equation (15) in the inhomogeneous matrix-inversion representation

$$\left({\mathcal{A}}^{-1}\left(\epsilon \right)+\lambda V^{\left(K\right)}{\mathsf{\Pi }}^{(M⨁N)}\right)\left(\begin{array}{c}|{\overrightarrow{y}}^{\left(a\right)}\succ \\ |{\overrightarrow{y}}^{\left(b\right)}\succ \end{array}\right)=\left(\begin{array}{c}|r^{\left(a\right)}\succ \\ |r^{\left(b\right)}\succ \end{array}\right)$$

(24)

or, equivalently,

$$\left(I+\lambda \mathcal{A}\left(\epsilon \right)V^{\left(K\right)}{\mathsf{\Pi }}^{(M⨁N)}\right)|\overrightarrow{y}\rangle =\mathcal{A}\left(\epsilon \right)|r\rangle .$$

(25)

Once we drop the redundant superscripts, and once we add the relevant parameter-dependences in (25) we obtain relation

$$\left(I+\lambda \mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\right)|\overrightarrow{y}\rangle =\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle .$$

(26)

This is our ultimate, iteration-friendly exact form of our perturbed Schrödinger equation.

4.3. Solutions

Equation (26) yields the ket-vector part of the solution in closed form,

$$|\overrightarrow{y}\rangle ={\left(I+\lambda \mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\right)}^{-1}\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle =|{\overrightarrow{y}}^{\left(solution\right)}(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle .$$

(27)

In the small-perturbation regime the latter formula may be given the conventional Taylor-series form with

$$\begin{gathered} |{\overrightarrow{y}}^{\left(solution\right)}(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle =\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle -\lambda \mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle + \\ +{\lambda }^2\mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle -\dots . \end{gathered}$$

(28)

Naturally, the construction is not yet finished because what is missing is the guarantee of equivalence between the eigenvalue problem (15) and its matrix-inversion reformulation (24) containing two redundant parameters. This is the last obstacle, easily circumvented by our setting, in solution (28), both of the redundant upper-case constants ${\mathsf{\Omega }}_{\left(a\right)}=y_M^{\left(a\right)}$ and ${\mathsf{\Omega }}_{\left(b\right)}=y_N^{\left(b\right)}$ equal to zero. In the light of explicit formula (28), this requirement is equivalent to the pair of relations

$$y_M^{\left(solution\right)}(\lambda ,\epsilon ,{\omega }_{\left(a\right)},{\omega }_{\left(b\right)})=0,y_{M+N}^{\left(solution\right)}(\lambda ,\epsilon ,{\omega }_{\left(a\right)},{\omega }_{\left(b\right)})=0.$$

(29)

Both left-hand-side functions of the three unknown quantities $\epsilon$, ${\omega }_{\left(a\right)}$ and ${\omega }_{\left(b\right)}$ are available, due to formula (28), in closed form. Although both can vary with the three independent unknowns (i.e., with $\epsilon ,{\omega }_{\left(a\right)}$ and ${\omega }_{\left(b\right)}$), one of these variables merely plays the role of an optional normalization constant so that we may set, say, ${\omega }_{\left(a\right)}^2+{\omega }_{\left(b\right)}^2=1$. Thus, the implicit version of the perturbation-expansion construction of bound states is completed.

5. Schrödinger Equation in Leading-Order Approximation

Two coupled equations (29) determine the bound state. In the spirit of perturbation theory one may expect that the perturbations happen to be, in some sense, small. At the same time, even the analysis of the comparatively elementary $L=1$ secular Equation (A16) determining the single free variable (viz., the energy) led to the necessity of a strongly counterintuitive scaling (A18) reflecting, near the extreme EP boundary of unitarity, the strong anisotropy of the geometry of the physical Hilbert space. Naturally, at least comparable complications must be expected to be encountered during the analysis of the more complicated set of two coupled equations (29) representing the secular equation in its exact $L=2$ version, constructed as particularly suitable for systematic approximations.

5.1. Generic Case: Perturbations without Vanishing Elements

Even the most drastic truncation of the formal power series (28) yields already a nontrivial ket vector

$$|{\overrightarrow{y}}^{\left(solution\right)}(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle =\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle .$$

(30)

Needless to add that what must vanish are the auxiliary variables ${\mathsf{\Omega }}_{(a,b)}$ alias two functions which are available in closed form. Thus, we must solve the following two simplified equations

$$\sum _{k=1}^{N_{\left(\varrho \right)}}{A}_{N_{\left(\varrho \right)},k}^{\left(\varrho \right)}\left(\epsilon \right)r_k^{\left(\varrho \right)}(\lambda ,\epsilon ,{\omega }_{\left(a\right)},{\omega }_{\left(b\right)})=0,\varrho =a,b.$$

(31)

After the insertion of the respective matrix elements $A_{N_{(a,b)},k}^{(a,b)}\left(\epsilon \right)$ [cf. Equation (A9)] we obtain the pair of relations

$$\left[\sum _{m=0}^{M-1}{ϵ}^m\lambda V_{M-m,1}^{(M,M)}\right]{\omega }_{\left(a\right)}+\left[\sum _{m=0}^{M-1}{ϵ}^m\lambda V_{M-m,1}^{(M,N)}\right]{\omega }_{\left(b\right)}={ϵ}^M{\omega }_{\left(a\right)},$$

(32)

$$\left[\sum _{n=0}^{N-1}{ϵ}^n\lambda V_{N-n,1}^{(N,M)}\right]{\omega }_{\left(a\right)}+\left[\sum _{n=0}^{N-1}{ϵ}^n\lambda V_{N-n,1}^{(N,N)}\right]{\omega }_{\left(b\right)}={ϵ}^N{\omega }_{\left(b\right)}.$$

(33)

Obviously, this set can be read as a generalized eigenvalue problem which determines generalized eigenvectors $\overrightarrow{\omega }$ at a K-plet of eigenenergies $ϵ$ which are all defined as roots of the corresponding generalized secular determinant.

5.2. Hierarchy of Relevance and Reduced Approximations

In the generic case one must assume that the matrix elements of the perturbation do not vanish and that $\lambda$ is small, i.e., that also the eigenvalues $\epsilon$ remain small. This enables one to omit all the asymptotically subdominant corrections and to consider just the linear algebraic system

$$\left[\lambda V_{M,1}^{(M,M)}-{ϵ}^M\right]{\omega }_{\left(a\right)}+\lambda V_{M,1}^{(M,N)}{\omega }_{\left(b\right)}=0,$$

(34)

$$\lambda V_{N,1}^{(N,M)}{\omega }_{\left(a\right)}+\left[\lambda V_{N,1}^{(N,N)}-{ϵ}^N\right]{\omega }_{\left(b\right)}=0.$$

(35)

The solution of these two simplified coupled linear relations exists if and only if the determinant of the system vanishes,

$$det\left[\begin{array}{cc}\lambda V_{M,1}^{(M,M)}-{ϵ}^M& \lambda V_{M,1}^{(M,N)}\\ \lambda V_{N,1}^{(N,M)}& \lambda V_{N,1}^{\left(NN\right)}-{ϵ}^N\end{array}\right]=0.$$

(36)

Thus, an ordinary eigenvalue problem is encountered when $M=N$.

Lemma 1.

In the generic equipartitioned cases with $M=N\ge 3$ the spectrum ceases to be all real under bounded perturbations. The loss of unitarity is encountered.

Proof. 

Both roots ${ϵ}^N=\lambda x$ of the exactly solvable quadratic algebraic secular Equation (36) may be guaranteed to be real in a certain domain of parameters. Still, some of the energies themselves are necessarily complex since $ϵ=\sqrt[N]{\lambda x}$ is an N-valued function with values lying on a complex circle.  □

In the other, non-equipartitioned dynamical scenarios with, say, $M>N$, the behavior of the system in an immediate vicinity of its EP extreme is still determined by the asymptotically dominant part of the secular equation. After an appropriate modification, the above proof still applies.

Lemma 2.

In the generic case with $M>N\ge 3$ we get, from the dominant part of the generalized eigenvalue problem (36), a subset (N-plet) of asymptotically dominant eigenvalues $\epsilon =\mathcal{O}({\lambda }^{1/N})$ which cannot be all real.

5.3. Unitary Case: Re-Scaled Perturbations

The loss of unitarity occurring in the $L=1$ models was associated with the use of the too broad a class of norm-bounded perturbations. From our preliminary results described in preceding subsection one can conclude that a similar loss of unitarity may be also expected to occur, in the generic case, at $L=2$. Indeed, the anisotropy of the physical Hilbert space which reflects the influence of an EP-related singularity of the Hamiltonian may be expected to lead again to a selective enhancement of the weight of certain specific matrix elements of perturbations $\lambda V^{\left(K\right)}$.

Once we recall the $L=1$ scenario of Appendix B.3 we immediately imagine that the main source of the apparent universality of the instability under norm-bounded perturbations should be sought, paradoxically, in the routine but, in our case, entirely inadequate norm-boundedness assumption itself. Indeed, in a way documented by Lemmas 1 and 2, the loss of the reality of spectra may directly be attributed to the conventional and comfortable but entirely random, unfounded and formal assumption of the uniform boundedness of the matrix elements of the perturbations.

Most easily the latter result may be illustrated using the drastically simplified version (36) of the leading-order secular equation. Dominant role is played there by the quadruplet of matrix elements $V_{(P,1)}^{(P,Q)}$ with superscripts P and Q equal to M or N. Once they are assumed $\lambda$-independent and non-vanishing, the leading-order energies read $ϵ=\sqrt[N]{\lambda x}$. Thus, at any $N\ge 3$ their N-plet forms an equilateral N-angle in the complex plane of $\lambda$.

The latter observation inspires a remedy. In a way eliminating the $N\ge 3$ complex-circle obstruction one simply has to re-scale the energies as well as all the relevant matrix elements of the perturbation. In this manner the ansatz

$$ϵ=ϵ\left(E\right)=\sqrt{\lambda }E$$

(37)

opens the possibility of the spectrum being real. Another multiplet of postulates

$$V_{M-m,1}^{(M,M)}={\lambda }^{(M-m)/2}{W}_{M-m,1}^{(M,M)},V_{M-m,1}^{(M,N)}={\lambda }^{(M-m)/2}{W}_{M-m,1}^{(M,N)},m=0,1,\dots ,M-1,$$

(38)

and

$$V_{N-n,1}^{(N,M)}={\lambda }^{(N-n)/2}{W}_{N-n,1}^{(N,M)},V_{N-n,1}^{(N,N)}={\lambda }^{(N-n)/2}{W}_{N-n,1}^{(N,N)},n=0,1,\dots ,N-1$$

(39)

now contains a new partitioned matrix $W=W^{\left(K\right)}$ which is assumed uniformly bounded.

Lemma 3.

There always exists a non-empty $(2M+2N)$-dimensional domain $\mathcal{D}$ of the “physical” matrix elements of W for which the leading-order spectrum is all real and non-degenerate, i.e., in the language of physics, tractable as stable bound-state energies.

Proof. 

The main consequence of the amended, necessary-condition perturbation-smallness requirements (38) and (39) is that in our initial, unreduced leading-order generalized eigenvalue problem (32) + (33), all terms in the sums become of the same order of magnitude. By the scaling we managed to eliminate the explicit presence of the measure of smallness $\lambda$. In other words the input information about dynamics is formed now by the re-scaled perturbation matrix W which offers an $(2M+2N)$-plet of free $\mathcal{O}\left(1\right)$ parameters. At the same time, the spectral-definition output is given by the $M+N$ roots of the corresponding secular equation, i.e., by the roots of an $(M+N)$-th-degree polynomial in E, with all its separate coefficients bounded and, in general, non-vanishing. Under these conditions the assertion of the lemma is obvious.  □

6. Discussion

6.1. Schrödinger Picture and Quasi-Hermitian Hamiltonians

In the context of the new, Bender- and Boettcher-inspired paradigm the non-Hermiticity of $H\left(\gamma \right)$ in $\mathcal{K}$ may be considered compatible with unitarity whenever the spectrum itself is found real. The explanation of the apparent paradox is easy: Under certain reasonable mathematical assumptions (cf. [21]) one can find a non-unitary invertible mapping $\mathsf{\Omega }=\mathsf{\Omega }\left(\gamma \right)$ with property ${\mathsf{\Omega }}^{†}\mathsf{\Omega }=\mathsf{\Theta }\ne I$ which makes our Hamiltonian self-adjoint,

$$H\left(\gamma \right)\to \mathfrak{h}\left(\gamma \right)=\mathsf{\Omega }\left(\gamma \right)H\left(\gamma \right)\mathsf{\Omega }\left(\gamma \right)={\mathfrak{h}}^{†}\left(\gamma \right).$$

(40)

In practice the idea enables one to avoid the use of the “conventional” Hamiltonian $\mathfrak{h}\left(\gamma \right)$ (self-adjoint, by construction, in another Hilbert space $\mathcal{L}$) whenever it happens to be “prohibitively complicated”. In the literature such a type of simplification of calculations is usually attributed to Dyson [22]. Expectedly, the strategy (also known as preconditioning) proved efficient as a tool of construction of bound states in nuclear physics [21]. In spite of a certain initial doubts [23,24], the approach also proved applicable in the quantum physics of scattering [25,26].

From a more abstract theoretical point of view the reference to the “Dyson’s” isospectral mapping (40) becomes redundant when we reclassify the space $\mathcal{K}$ as “unphysical”, and when we redefine its inner product yielding another, “physical” Hilbert space $\mathcal{H}$,

$${\langle {\psi }_1|{\psi }_2\rangle }_{\mathcal{H}}=\langle {\psi }_1\left|\mathsf{\Theta }\right|{\psi }_2{\rangle }_{\mathcal{K}},|{\psi }_{1,2}\rangle \in \mathcal{K}.$$

(41)

A new, equivalent, “two-Hilbert-space” version of the Schrödinger picture is obtained in which the physics described in the correct Hilbert space $\mathcal{H}$ is “translated” to its mathematically easier representation in $\mathcal{K}$. In this sense, the self-adjointness of hypothetical $\mathfrak{h}$ in hypothetical $\mathcal{L}$ is found equivalent to the self-adjointness of H in $\mathcal{H}$, represented by the relation

$$H^{†}\left(\gamma \right)\mathsf{\Theta }\left(\gamma \right)=\mathsf{\Theta }\left(\gamma \right)H\left(\gamma \right).$$

(42)

Dieudonné [27] and Scholtz et al [21] suggested to call relation (42) the quasi-Hermiticity of H in $\mathcal{K}$. In the context of quantum phenomenology a decisive amendment of the formalism may be seen in the split of the description of dynamics with information carried by the two operators $H\left(\gamma \right)$ and $\mathsf{\Theta }\left(\gamma \right)$ defined in $\mathcal{K}$ in place of one (viz., of $\mathfrak{h}\left(\gamma \right)$ living in $\mathcal{L}$).

6.2. Non-Hermitian Degeneracies with $L>2$

After the present decisive clarification of the possibility of the replacement of the $L=1$ formalism by its $L=2$ generalization, the next move to the further, $L>2$ dynamical scenarios is now an almost elementary exercise. Indeed, in Section 3, Section 4 and Section 5 it would be sufficient to move from the $L=2$ partitioning of $K=N_{\left(a\right)}+N_{\left(b\right)}$ to its arbitrary $L>2$ analogues (A1) and to the related partitioned vector sets [sampled, e.g., by Equation (7)] and matrices [sampled, e.g., by Equation (8)]. The $L=2$ doublets of the eligible superscripts ${}^{\left(a\right)}$ and ${}^{\left(b\right)}$ marking the curly kets [cf. (17) or (31), etc.] may be very easily extended to the L-plets with $L>2$, etc.

Along these lines the form of our basic power-series expansion of the perturbed bound-state kets (28) remains unchanged. The related $L=2$ compatibility constraint (29) must only be replaced by an L-plet of equations

$$y_{N_{(a_j)}}^{\left(solution\right)}(\lambda ,\epsilon ,\overrightarrow{\omega })=0,j=1,2,\dots ,L$$

(43)

where the unknown vector $\overrightarrow{\omega }$ has L components.

6.3. The Next-to-Leading-Order Approximation

In our present paper we did not pay too much attention to the next-to-leading-order (NLO) approximation where one would have to set

$$|{\overrightarrow{y}}^{\left(solution\right)}(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle =\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle -\lambda \mathcal{A}\left(\epsilon \right)V\mathsf{\Pi }\mathcal{A}\left(\epsilon \right)|r(\lambda ,\epsilon ,\overrightarrow{\omega })\rangle .$$

(44)

Our main reason was that such a generalization would make the associated compatibility conditions (43) perceivably more complicated. Indeed, in contrast to the leading-order case, an almost inessential improvement of the insight in the qualitative features of the quantum system in question would be accompanied, in the NLO formulae, by the emergence of multiple new matrix elements of perturbations $\lambda V^{\left(K\right)}$ including even the terms which would be quadratic functions of these matrix elements.

This being said, it is necessary to add that even on the pragmatic and qualitative level, the omission of the NLO corrections only seems completely harmful in the open-system setting using bounded matrices $\lambda V^{\left(K\right)}$ where one does not insist on having the strictly real spectrum. The point is that whenever one must guarantee the unitarity of the system, the class of the admissible perturbations must be further restricted. In this sense, unfortunately, an analysis using NLO might prove necessary. Indeed, the use and precision of the leading-order approximation need not be sufficiently reliable in general. Even a quick glimpse at the underlying assumptions (38) and (39) reveals that the leading-order approximation does not incorporate the influence of a large subset of the matrix elements of $\lambda V^{\left(K\right)}$. At the same time, one must keep in mind that in our present approach the dimension of the matrices K has been assumed finite. For this reason, in the case of doubts, a turn to the more universal and brute-force numerical methods might prove to be, in practical calculations, a reasonable alternative to the rather lengthy and complicated NLO calculations.