Adiabatic state conversion for (a)cyclic non-Hermitian quantum Hamiltonians of generalized functional form

Nye, Nicholas S.; Kantartzis, Nikolaos V.

doi:10.1063/5.0225403

In the present paper, we aim to firmly establish the adiabatic properties of two-level non-Hermitian quantum structures evolving along generalized (open/acyclic or closed/cyclic) paths in parameter space. Analytical solutions in terms of Airy and modified Bessel functions have been retrieved for linear and hyperbolic temporal dependencies in parity-time-symmetric-like systems, which were subsequently studied in the slowly varying limit to show conversion to one of the instantaneous eigenstates. Such a mode switching behavior is found to be an identifying feature of dissipative quantum settings, whether they evolve along cyclic or acyclic trajectories, and this has been proven in our paper by separately analyzing the dynamics of (i) the ratio of the state vector components, via a variant of the Möbius transformation, and (ii) the complex probability amplitudes, through a systematic inspection of the mode population equations. In the latter instance, it was furthermore shown that the identity of the eigenstate, to which the quantum arrangement transitions, depends highly on the magnitude of the adiabatic rate β. Along these lines, the concepts of the instantaneous $(D)$ and averagely $(D_{a v})$ dominant eigenstates are brought forth, while a reconfigurable photonic switch is also proposed, which can convert either to the $D$ or to the $D_{a v}$ modes based on the total period of evolution. Finally, we apply our findings in the case of closed parametric paths to demystify the recently reported symmetric and asymmetric state conversion effects and additionally demonstrate that operation at or near exceptional points does not qualitatively affect the conclusions of the current investigation.

I. INTRODUCTION

Non-Hermitian quantum mechanics constitutes an important extension to the standard Hermitian Hamiltonian formalism, as it enables us to investigate the dynamics of systems (S) that interact with the environment (E) and thus exhibit quantum dissipation. In this case, we say that the arrangement under study (S) suffers from decoherence (i.e., a definite phase relation no longer exists between the different quantum states), which in turn leads to information loss (from S to E) and to the generation of entanglement between the system and its surroundings. To overcome this issue and obtain a problem that can be still described by Hermitian Hamiltonians, we could have simply considered S and E together as a single unit. Following this thought process incrementally, it would eventually require having knowledge of the state of the entire universe, which even desirable is not possible. For practical purposes, instead of focusing on an infinite number of degrees of freedom, it is preferable to limit our attention to a smaller number of relevant variables describing the subsystem S and subsequently attain the corresponding effective equations of motion. The latter take into account any irrelevant degrees of freedom (associated with the interaction between S and E) via appropriate stochastic or dissipative terms. It should also be noted that decoherence can arise not only due to energy exchange with the environment but as a direct consequence of the act of measuring S (quantum decoherence/measurement theory/dissipation/entanglement are all highly interrelated concepts^1,2). Asides from all the previous discussion on the importance of open quantum arrangements, it is noteworthy to mention that certain quantum dynamical problems (e.g., decay or capture resonance effects and electron/photon scattering from atoms or molecules) can be more efficiently and comprehensively treated when employing outgoing/incoming wave boundary conditions, which by itself demands the definition of non-Hermitian Hamiltonians³ (their Hermitian counterparts require that the wavefunction vanishes at the boundary/infinity). Of course, the applicability of non-Hermitian theory is by no means limited to the quantum domain but rather extended to a plethora of other fields dealing with non-conservative or stochastic settings, including biological physics, photonics, optomechanics, hydrodynamics, and neural networks.^3,4

In recent years, a flurry activity has been noticed in the field of non-Hermitian physics and this has been greatly incited from the seminal paper by Bender and Boettcher⁵ on parity-time $(P T)$ -symmetric Hamiltonians (i.e., Hamiltonians H that commute with the parity-time operator according to $[P T, H] = 0$ ⁠). Such non-conservative settings exhibit a spontaneous symmetry breaking point, below and above which the eigenspectra become entirely real (the $P T$ symmetry of Hamiltonian H is then said to be unbroken, given that $[P T, H] = 0$ guarantees that $P T$ and H share a simultaneous set of eigenstates) and entirely complex (the $P T$ symmetry of H is then broken, given that $[P T, H] = 0$ no longer implies that $P T$ and H share the same eigenstates owing to the antilinearity of $T$ —parity operator $P$ ⁠: linear; time reversal operator $T$ ⁠: antilinear), respectively. An even greater peculiarity is that at this $P T$ symmetry breaking threshold, eigenvectors and eigenvalues coalesce; such a degeneracy, which is a characteristic attribute of non-Hermitian structures, is known in the literature as an exceptional point (EP) and has led to a number of intriguing observations.^6–9 Latest research has focused on the slowly varying properties of these open Hamiltonian arrangements. According to the standard quantum adiabatic theorem, as originally developed by Born and Fock^10–12 in the context of Hermitian theory, a physical setting is expected to remain in its instantaneous eigenstate assuming that any external influences act on it slowly enough. Yet the previous statement is violated within non-conservative environments, where it has been shown that mode switching occurs for evolution along cyclic paths in the Hamiltonian parameter space^13–22 and this has been exploited for the demonstration of asymmetric transport in a variety of classical and quantum experimental setups.^23–28 Here, we extend these findings for acyclic parametric trajectories and more importantly unveil a previously overlooked fact: the identity of the mode describing the resulting quantum state highly depends on the magnitude of the adiabatic rate of evolution β. Capitalizing on the latter, we shall elaborate on the results reported in cyclic non-Hermitian configurations and also clear out any prior misconceptions regarding the correlation between EP encirclement and symmetric (SSC)/asymmetric (ASC) effects—more specifically, it will be seen that encircling EPs is neither a necessary nor a sufficient condition for SSC or ASC to take place.

Overall, the present article establishes the fundamental principles of open quantum systems, pertaining to their adiabatic evolution behavior (in what follows, the terms adiabatic, slowly varying, and quasistatic shall be used equivalently and interchangeably). While here we emphasize on two-level Hamiltonians, the physical insights and mathematical reasoning can be extended to higher-dimensional Hilbert spaces (particular examples related to 3 × 3 non-Hermitian Hamiltonians are included in the supplementary material; however, this topic will be explored in greater depth in our future research). Along these lines, the dynamical equations dictating the quantum state vector evolution associated with a 2 × 2 generalized and time-dependent Hamiltonian are initially retrieved (Sec. II A). Appropriate state vector transformations are also discussed, with the purpose of broadening the range of analytically solvable scenarios. Via the transfer matrix method, it is shown how the full system dynamics can be attained at any time instant, and this is subsequently applied in the special case of $P T$ -symmetric-like Hamiltonian arrangements (Sec. II B). The parameter space is then defined by the gain/loss and detuning factors, and in this respect, two different types of parametric paths were examined: linear (this leads to parabolic cylinder and Airy solutions, with the latter applicable under certain approximations—Sec. II B 1) and hyperbolic (the dynamics can then be expressed in terms of modified Bessel functions—Sec. II B 2). The slowly varying response was studied in both circumstances, by employing the asymptotic properties of the corresponding solutions. Under an acyclic variation of the system parameters, it was illustrated that the resulting quantum state is always described by one of the instantaneous eigenstates of the non-Hermitian Hamiltonian, regardless of any initial excitation conditions. This phenomenon, known as state conversion or mode switching, is intimately linked to the Stokes phenomenon of asymptotics, as will also become evident through our detailed mathematical calculations.

In Secs. III and IV, we are concerned with the quasistatic features of generalized non-Hermitian settings [the respective trajectories within the parameter space can assume arbitrary geometric forms and may be either open (Sec. III) or closed (Sec. IV)]. In this regard, our study considers separately the evolution of (i) the ratio of the state vector components (Sec. III) and (ii) the mode population amplitudes (Sec. III and Appendix B 3). A discretization of the Schrödinger equation into small time steps indicates that quantity (i) obeys a time-dependent version of the Möbius transformation. By leveraging on the properties of the latter (after an appropriate adaptation for the time-dependent case), it is shown that under non-Hermitian and sufficiently slowly varying conditions, the quantum state converts to the least dissipative eigenstate except for very special circumstances (in the Hermitian case, the standard quantum adiabatic theorem is confirmed). Our derivation generalizes the findings of our recent work,²² where it is theoretically demonstrated that for small time steps, state conversion occurs in discrete non-Hermitian structures for a cyclic variation of the system parameters.

Regarding now quantity (ii), the quantum state is expressed at each time instant in terms of the instantaneous eigenmodes, in a manner also consistent with Wentzel–Kramers–Brillouin (WKB) theory for very large evolution periods (the WKB theory results were additionally utilized to confirm the asymptotic behavior of the analytical solutions retrieved in Secs. II B 1 and II B 2, in the context of $P T$ -symmetric-like configurations). A rigorous investigation of the mode population equations then showed that the resulting system dynamics will be governed by a single eigenstate, whose identity depends on the degree of adiabaticity (the role of the non-adiabatic mode couplings has also been clearly highlighted). In this respect, the concepts of instantaneous/averagely dominant eigenstates were employed to describe moderately/extremely slow non-Hermitian processes with a special emphasis on the topology of the associated eigenspectra evolution graphs. This observation was subsequently used to propose a reconfigurable optical switch based on discrete photonic lattices, which can operate in either a cyclic or an acyclic fashion—in such a classical optics platform, the state is defined through the polarization of the propagating field, its dynamics are dictated by a discrete Schrödinger-like equation, while the adiabatic rate of evolution is controlled via the number of round trips along an appropriately designed fiber loop (for more technical details on the Schrödinger-type equations governing discrete or continuous photonic setups, see Refs. 22, 28, and 29). The case of closed paths in parameter space was examined in a separate section (Sec. IV) in order to (a) provide a more detailed interpretation of the symmetric/asymmetric mode switching effects observed within non-conservative environments^13–28 (here, in essence, we apply the analytic findings of Sec. III pertaining to parametric trajectories of arbitrary geometric form) and (b) demonstrate the functional characteristics of the aforementioned optical switch, showcasing its capacity to support a dual operational regime (SSC or ASC) determined by the value of β. Finally, it has been shown that the crossing of spectrum degeneracies can have only a quantitative (magnitude of β for conversion to take place) rather than a qualitative impact in our study.

II. ASYMPTOTIC ANALYSIS OF EXACTLY SOLVABLE OPEN QUANTUM SYSTEMS

A. Two-level Hamiltonian settings

The dynamics of a generalized linear bipartite Hamiltonian system can be described in terms of the continuous time Schrödinger equation as

i \frac{\partial |u (t)⟩}{\partial t} = H (t) |u (t)⟩ = [\begin{matrix} h^{11} (t) & h^{12} (t) \\ h^{21} (t) & h^{22} (t) \end{matrix}] |u (t)⟩,

(1)

where t ∈ [0, T] (T is the overall evolution time) reflects the temporal variable (the above expression is exact for optical setups, but for quantum mechanical arrangements, the substitution t → t/ℏ should be performed to be consistent with the units of the various quantities), H(t) is the time-dependent quantum Hamiltonian, and

|u (t)⟩ = {[a (t) b (t)]}^{T}

(⁠

T

indicates the matrix transposition operator) represents the state vector. Here, we shall use the following equivalent form to Eq. :

\frac{\partial |u (t)⟩}{\partial t} = M (t) |u (t)⟩ = [\begin{matrix} m^{11} (t) & m^{12} (t) \\ m^{21} (t) & m^{22} (t) \end{matrix}] |u (t)⟩,

(2)

which is commonly employed in dynamical systems theory and also eliminates the need to be concerned about the additional factor i (⁠

i = \sqrt{- 1}

stands for the imaginary unit). To distinguish between M(t) and H(t), we shall call the former as the system matrix and it will evidently apply that M(t) = −iH(t). Given that in the present work, we are particularly concerned with examining the quasistatic properties of Eqs. and , it will be assumed that matrices H(t) and M(t) can be written as functions of the scaled time coordinate

\bar{t} = t / T

[i.e.,

H = H (\bar{t}), M = M (\bar{t})

] with β = 1/T defining the adiabatic rate. It will be true then for adiabatic evolution conditions (T → ∞) that

d Q (t) / d t = T^{- 1} d Q (\bar{t}) / d \bar{t} ≃ 0

⁠, where Q depicts the matrices H, M, or their respective eigenvalues/eigenvectors.

Conventional Hermitian arrangements satisfy the condition H^†(t) = H(t) or equivalently M^†(t) = −M(t) [Q^†(t) refers to the conjugate transpose/adjoint of Q(t)], which in turn demand that the eigenvalues/eigenenergies $λ_{H}^{\pm} (t)$ of the Hamiltonian matrix H(t) are purely real and that the eigenvalues λ^±(t) $[= - i λ_{H}^{\pm} (t)]$ of the system matrix M(t) are purely imaginary [symbols (+) and (−) designate the modes of the two-level Hamiltonian structure under study]. Here, we are particularly interested in the non-Hermitian case (open quantum systems), where the eigenvalues $λ_{H}^{\pm} (t), λ^{\pm} (t)$ are in general complex. The imaginary (real) part of $λ_{H}^{\pm} (t)$ [λ^±(t)] may signify (i) gain/loss in optical setups and (ii) dissipation/decoherence in quantum mechanical setups. Existing analytical or experimental investigations into state conversion effects in slowly varying non-conservative settings^13–28 predominantly rely on a cyclic variation of the system parameters. In the current paper, such a condition is relaxed, i.e., we allow the Hamiltonian arrangement to move along open paths $C^{R}$ in parameter space $R$ ⁠.

Since throughout the analysis, an explicit temporal dependence of the various Hamiltonian parameters is assumed, the partial derivative (∂) in Eqs. and can be replaced with a total derivative (d). As such, the linear differential equation that the state vector component a(t) satisfies can be found to be [Eq. ]

\begin{align} \frac{d^{2} a (t)}{d t^{2}} - \{t r [M (t)] + \frac{d m^{12} (t)}{d t} / m^{12} (t)\} \frac{d a (t)}{d t} \\ + \{\det [M (t)] + m^{11} (t) \frac{d m^{12} (t)}{d t} / m^{12} (t) - \frac{d m^{11} (t)}{d t}\} a (t) = 0, \end{align}

(3)

where tr[M(t)] = m¹¹(t) + m²²(t) and det[M(t)] = m¹¹(t)m²²(t) − m¹²(t)m²¹(t) denote the trace and determinant of the time-dependent matrix M(t), respectively [the evolution equation for b(t) can be obtained after a simple exchange of superscripts 1 ↔ 2 in the above expression]. The solution to Eq. can be written as a linear superposition of two linearly independent functions in the form

a (t) = c_{1}^{a} a_{1} (t) + c_{2}^{a} a_{2} (t)

⁠, with

c_{1, 2}^{a}

indicating constants and the Wronskian of a₁(t), a₂(t) given as (Abel’s identity)

\begin{align} W (t) & = a_{1} (t) \frac{d a_{2} (t)}{d t} - a_{2} (t) \frac{d a_{1} (t)}{d t} \\ = C^{W} \frac{m^{12} (t)}{m^{12} (0)} e^{\int_{0}^{t} t r [M (t^{'})]} d t^{'} \neq 0 . \end{align}

(4)

The presence of constant

C^{W}

[=W(0) ≠ 0] in the previous formula allows for the arbitrary choice of a lower time limit in the integral appearing in the exponent. The combined dynamics of a(t), da(t)/dt can be expressed in the following matrix form:

[\begin{matrix} a (t) \\ \frac{d a (t)}{d t} \end{matrix}] = Z (t) [\begin{matrix} c_{1}^{a} \\ c_{2}^{a} \end{matrix}] = [\begin{matrix} a_{1} (t) & a_{2} (t) \\ \frac{d a_{1} (t)}{d t} & \frac{d a_{2} (t)}{d t} \end{matrix}] [\begin{matrix} c_{1}^{a} \\ c_{2}^{a} \end{matrix}] .

(5)

The constants

c_{1, 2}^{a}

can be directly determined from the initial conditions for a(t) and da(t)/dt, and as such, it will apply that

[\begin{matrix} a (t) \\ \frac{d a (t)}{d t} \end{matrix}] = Z (t) Z^{- 1} (0) {[\begin{matrix} a (t) \\ \frac{d a (t)}{d t} \end{matrix}]}_{0} = Z (t, 0) {[\begin{matrix} a (t) \\ \frac{d a (t)}{d t} \end{matrix}]}_{0},

(6)

where the subscript 0 indicates the evaluation at t = 0 {matrix

Z (t)

is invertible since

\det [Z (t)] = W (t) \neq 0

}. The elements of transfer matrix Z(t, 0) will be given by

ζ^{11} (t, 0) = \frac{a_{1} (t) {[\frac{d a_{2} (t)}{d t}]}_{0} - a_{2} (t) {[\frac{d a_{1} (t)}{d t}]}_{0}}{W (0)},

(7a)

ζ^{12} (t, 0) = \frac{a_{2} (t) a_{1} (0) - a_{1} (t) a_{2} (0)}{W (0)},

(7b)

ζ^{21} (t, 0) = \frac{\frac{d a_{1} (t)}{d t} {[\frac{d a_{2} (t)}{d t}]}_{0} - \frac{d a_{2} (t)}{d t} {[\frac{d a_{1} (t)}{d t}]}_{0}}{W (0)},

(7c)

ζ^{22} (t, 0) = \frac{\frac{d a_{2} (t)}{d t} a_{1} (0) - \frac{d a_{1} (t)}{d t} a_{2} (0)}{W (0)} .

(7d)

We also state the following useful formula, which relates the ratio of the state vector components [r(t) = b(t)/a(t)] to [da(t)/dt]/a(t) [see Eq. ]:

r (t) = \frac{1}{m^{12} (t)} \frac{d a (t) / d t}{a (t)} - \frac{m^{11} (t)}{m^{12} (t)} .

(8)

Equations – allow us now to evaluate the dynamics associated with a(t), da(t)/dt, and r(t) [and thus also b(t)], given a set of initial conditions. It should also be highlighted that the dynamical behavior of r(t) can be retrieved directly from Eq. as

\frac{d r (t)}{d t} + m^{12} (t) r^{2} (t) + [m^{11} (t) - m^{22} (t)] r (t) - m^{21} (t) = 0,

(9)

which corresponds to a nonlinear Riccati differential equation. Depending on the problem at hand and the resulting form of the linear [Eq. ] and nonlinear differential relations [Eq. ], either Eqs. – or Eq. can be employed to predict the temporal evolution of r(t).

Finally, we shall introduce an alternative formulation of the two-mode Hamiltonian system under study, which will enable us to identify a wider range of parametric temporal dependencies leading to analytically solvable scenarios. In this regard, the subsequent transformation will be applied on the state vector components,

\tilde{a} (t) = e^{- \int_{0}^{t} m^{11} (t^{'}) d t^{'}} a (t), \tilde{b} (t) = e^{- \int_{0}^{t} m^{22} (t^{'}) d t^{'}} b (t),

(10)

which was judiciously selected in order to recast Eqs. and into a hollow matrix form [in other words, the transformed components

\tilde{a} (t), \tilde{b} (t)

will experience solely a cross-modulation rather than a self-modulation effect]. Indeed, in this new framework, the dynamics will be described by

d |\tilde{u} (t)⟩ / d t = \tilde{M} (t) |\tilde{u} (t)⟩

⁠, where

|\tilde{u} (t)⟩ = {[\tilde{a} (t) \tilde{b} (t)]}^{T}

and the matrix

\tilde{M} (t) = - i \tilde{H} (t)

will be given by

\tilde{M} (t) = [\begin{matrix} 0 & m^{12} (t) e^{\int_{0}^{t} [m^{22} (t^{'}) - m^{11} (t^{'})] d t^{'}} \\ m^{21} (t) e^{\int_{0}^{t} [m^{11} (t^{'}) - m^{22} (t^{'})] d t^{'}} & 0 \end{matrix}] .

(11)

The choice of the lower limit in the integrals appearing in Eq. [and subsequently in Eq. ] implies that

\tilde{a} (0) = a (0)

and

\tilde{b} (0) = b (0)

⁠; however, any alternative selection that facilitates our analysis may be made. The linear differential dynamics governing

\tilde{a} (t)

can be found to be

\begin{align} \frac{d^{2} \tilde{a} (t)}{d t^{2}} - \{/ m^{22} (t) - m^{11} (t) + \frac{d m^{12} (t)}{d t} m^{12} (t)\} \frac{d \tilde{a} (t)}{d t} \\ - m^{12} (t) m^{21} (t) \tilde{a} (t) = 0, \end{align}

(12)

with the solution expressed as

\tilde{a} (t) = c_{1}^{\tilde{a}} {\tilde{a}}_{1} (t) + c_{2}^{\tilde{a}} {\tilde{a}}_{2} (t)

(⁠

c_{1, 2}^{\tilde{a}}

denote constants) and the Wronskian of

{\tilde{a}}_{1} (t), {\tilde{a}}_{2} (t)

given from Abel’s identity as [

C^{\tilde{W}} = \tilde{W} (0) \neq 0

depicts a constant in the succeeding relation]

\begin{align} \tilde{W} (t) & = {\tilde{a}}_{1} (t) \frac{d {\tilde{a}}_{2} (t)}{d t} - {\tilde{a}}_{2} (t) \frac{d {\tilde{a}}_{1} (t)}{d t} \\ = C^{\tilde{W}} \frac{m^{12} (t)}{m^{12} (0)} e^{\int_{0}^{t} [m^{22} (t^{'}) - m^{11} (t^{'})] d t^{'}} \neq 0 . \end{align}

(13)

The evolution of

\tilde{b} (t)

can be directly attained from Eq. via an exchange of indices 1 ↔ 2, and formulas analogous to Eqs. – can be similarly derived for the transformed system.

B. Traceless and symmetric non-Hermitian Hamiltonian model

Here, we shall assume a Hamiltonian matrix H(t) with components of the form h¹¹(t) = −h²²(t) = if(t) = ig(t) + σ(t) and h¹²(t) = h²¹(t) = −κ [m¹¹(t) = −m²²(t) = f(t) = g(t) − iσ(t), m¹²(t) = m²¹(t) = iκ], where g(t), σ(t), and κ are real-valued parameters [in a photonic setup, g(t) and σ(t) will represent the time-dependent gain/loss and detuning/phase mismatch factors, respectively, while κ will stand for the constant and (without loss of generality) positive coupling strength between the state vector components a(t) and b(t)]. The resulting matrix equation describing the underlying dynamics will be given by

\frac{d |u (t)⟩}{d t} = [\begin{matrix} f (t) & i κ \\ i κ & - f (t) \end{matrix}] |u (t)⟩,

(14)

which corresponds to an almost parity-time-symmetric arrangement [exact

P T

symmetry demands that σ(t) = 0 and g(t) = g(T_o − t) with T_o denoting the time reversal point in the definition of the time reversal operator—for a more extensive discussion, see the appendixes of Ref. 22]. The evolution of a(t) and

\tilde{a} (t) = e^{- \int_{0}^{t} f (t^{'}) d t^{'}} a (t)

is determined by

\frac{d^{2} a (t)}{d t^{2}} + [κ^{2} - f^{2} (t) - \frac{d f (t)}{d t}] a (t) = 0,

(15a)

\frac{d^{2} \tilde{a} (t)}{d t^{2}} + 2 f (t) \frac{d \tilde{a} (t)}{d t} + κ^{2} \tilde{a} (t) = 0,

(15b)

while analogous equations can be retrieved for b(t) and

\tilde{b} (t) = e^{\int_{0}^{t} f (t^{'}) d t^{'}} b (t)

after setting f(t) → −f(t) in the previous expressions. The ratio variable r(t) will evolve according to

\frac{d r (t)}{d t} + i κ r^{2} (t) + 2 f (t) r (t) - i κ = 0

(16)

and is related to [da(t)/dt]/dt as

r (t) = - i \frac{1}{κ} \frac{/ d a (t) d t}{a (t)} + i \frac{f (t)}{κ} .

(17)

It is interesting now to note that Eq. (15b) corresponds to a harmonic oscillator exhibiting a complex and time-dependent damping factor [2f(t)] with a positive or a negative real profile, depending on the sign of g(t) = Re[f(t)] (Re[·] and Im[·] will designate the real and imaginary parts of the associated complex quantities, respectively). This observation is essential, since it enables us to take advantage of the large amount of analytical studies (and known analytical solutions) for damped harmonic oscillator systems, by simply identifying the right temporal dependencies for f(t). In this regard, if we allow f(t) to vary linearly [f(t) = η_f,1t + η_f,0 with η_f,1, η_f,0 constants] or hyperbolically [ $f (t) = {(η_{f, 1} t + η_{f, 0})}^{- 1}$ ] with time, then Eq. (15b) can be immediately recognized to be of the Hermite- and Bessel-type, respectively [in Ref. 22, an even richer landscape of dynamics is attained for alternative temporal dependencies of f(t)]. In what follows, both of the aforementioned scenarios shall be considered in the asymptotic regime—more specifically, Eq. (15a) shall be used when f(t) = η_f,1t + η_f,0 (this will lead to the less familiar parabolic cylinder functions and under certain assumptions to the Airy functions) and Eq. (15b) shall be employed when $f (t) = {(η_{f, 1} t + η_{f, 0})}^{- 1}$ ⁠.

So far, our discussion has been purely theoretical, but it should be highlighted that time-dependent non-Hermitian models of the form shown in Eq. (14) have been realized in a variety of experimental settings.^23–28 What is even more noteworthy is the demonstrated capability of precisely tailoring the complex damping factor via electro-optic amplitude (tuning of the gain coefficient) and phase (tuning of the phase mismatch) modulator components within judiciously designed fiber loop networks²⁸ (in this case, the state of the system is described by the polarization of the propagating field). This in turn opens up the possibility of emulating a large range of physical arrangements obeying dynamics analogous to Eqs. (14)–(17), including the structures that will be analyzed in Secs. II B 1 and II B 2.

1. Linear time dependence for complex damping factor f(t)

We shall examine now the instance where both gain and detuning parameters depend linearly on time, as depicted in Figs. 1(a)–1(c). The factor f(t) = g(t) − iσ(t) can be expressed then as f(t) = η_f,1t + η_f,0, with parameters η_f,1 and η_f,0 evaluated as follows [the top and bottom signs will signify the paths AO and BO, respectively, with A (κ, ɛκ), B (κ, −ɛκ), O (κ − ɛκ, 0)]:

η_{f, 1} = \frac{- ε κ \pm i ε κ}{T} = - \frac{ε κ}{\sqrt{2} T} e^{\mp i π / 4},

(18a)

η_{f, 0} = κ \mp i ε κ ≃ κ e^{\mp i ε} .

(18b)

In the above formulas, we assumed small values of

ε (e^{i ε} ≃ 1 + i ε, ε \in R^{+})

⁠, which geometrically implies that path

C^{R}

will lie closer to the point (g, σ) = (κ, 0) (this is the exceptional point degeneracy characterizing the

P T

-symmetric-like system under study—see Appendix A 1). Moreover, to simplify our notation in both Secs. II B 1 and II B 2, the various time-dependent quantities q(t) will simply appear as q (the value of q at the start t = 0 and end t = T of the parametric evolution will be signified as q₀ and q_T, respectively), while in all our calculations, the principal root definition shall be employed {i.e., the phases of complex numbers will be provided in the interval (−π, π]—this aligns also with the definition of the imaginary unit as

i = \sqrt{- 1} = \sqrt{e^{i π}} = e^{i π / 2}

}.

FIG. 1.

View large Download slide

Slowly varying dynamics of the state vector component ratio r(t) for two distinct time dependencies of the complex damping factor: (a)–(e) f(t) = η_f,1t + η_f,0 (Sec. II B 1) and (f)–(j) $f (t) = {(η_{f, 1} t + η_{f, 0})}^{- 1}$ (Sec. II B 2). In the above (and the following) schematics, the gain and detuning coefficients are normalized with respect to the coupling constant $(\overset{̌}{g} = g / κ, \overset{̌}{σ} = σ / κ, \overset{̌}{f} = f / κ)$ ⁠, while the time with respect to the overall period of evolution $(\bar{t} = t / T, {\bar{t}}^{'} = t / T^{'})$ —see the last bullet points in Appendix A 1 for the rationale behind such normalizations. The first column [(a) and (f)] depicts the trajectory followed in parameter space [either AO (dotted lines) or BO (solid lines)—here, we chose ɛ = 0.5 and thus the beginning and ending values for $\overset{̌}{f}$ will be 1 ∓ iɛ = 1 ∓ i0.5 (top/bottom sign: AO/BO) and 1 − ɛ = 0.5, correspondingly], while the location of the EP degeneracy at $(\overset{̌}{g}, \overset{̌}{σ}) = (1, 0)$ is signified with a magenta asterisk. The second column [(b), (c), (g), and (h)] indicates the time dependency of variables $\overset{̌}{g} = R e (\overset{̌}{f})$ and $\overset{̌}{σ} = - I m (\overset{̌}{f})$ [in (b) and (g), the dotted and solid lines overlap as in both cases, the gain parameter varies with time in exactly the same manner for paths AO and BO], while the third [(d) and (i)] and fourth [(e) and (j)] columns depict the temporal dynamics of the ratio factor r for two different evolution periods [κT = 30 (dark colored curves in the third column), κT′ = 5 (light colored curves in the fourth column), and T′ = T/6, green and orange lines denote the real and imaginary elements of r, respectively]. The filled and empty triangles (circles) at $\bar{t} = {\bar{t}}^{'} = 1$ in the last two columns describe the real and imaginary parts, respectively, of the component ratio associated with eigenstate (−) [(+)]—it can be clearly seen now how switching either to eigenmode (−) or (+) occurs for paths AO and BO, correspondingly, for sufficiently slowly varying conditions (column 3), thus confirming the analysis of Secs. II B 1 and II B 2. Initial excitation was randomly set, but kept the same between panels (d) and (e) panels (i) and (j) for a given parametric trajectory (AO or BO), to definitively illustrate that state conversion is an artifact of the magnitude of the total evolution period, rather than some other cause. In all instances, dotted/solid and dark/light lines indicate results relevant to paths AO or BO and periods T or T′, respectively, while the coupling constant was fixed in all simulation scenarios. It should also be noted that ɛ = 0.5 might not exactly satisfy the condition ɛ ≪ 1 used in Secs. II B 1 and II B 2 (this assumption was made to simplify our asymptotic calculations and to more efficiently convey the mathematical mechanism behind the observed mode switching effect) but was selected as such to demonstrate that state conversion emerges independently of the scale of ɛ.

FIG. 1.

View large Download slide

Slowly varying dynamics of the state vector component ratio r(t) for two distinct time dependencies of the complex damping factor: (a)–(e) f(t) = η_f,1t + η_f,0 (Sec. II B 1) and (f)–(j) $f (t) = {(η_{f, 1} t + η_{f, 0})}^{- 1}$ (Sec. II B 2). In the above (and the following) schematics, the gain and detuning coefficients are normalized with respect to the coupling constant $(\overset{̌}{g} = g / κ, \overset{̌}{σ} = σ / κ, \overset{̌}{f} = f / κ)$ ⁠, while the time with respect to the overall period of evolution $(\bar{t} = t / T, {\bar{t}}^{'} = t / T^{'})$ —see the last bullet points in Appendix A 1 for the rationale behind such normalizations. The first column [(a) and (f)] depicts the trajectory followed in parameter space [either AO (dotted lines) or BO (solid lines)—here, we chose ɛ = 0.5 and thus the beginning and ending values for $\overset{̌}{f}$ will be 1 ∓ iɛ = 1 ∓ i0.5 (top/bottom sign: AO/BO) and 1 − ɛ = 0.5, correspondingly], while the location of the EP degeneracy at $(\overset{̌}{g}, \overset{̌}{σ}) = (1, 0)$ is signified with a magenta asterisk. The second column [(b), (c), (g), and (h)] indicates the time dependency of variables $\overset{̌}{g} = R e (\overset{̌}{f})$ and $\overset{̌}{σ} = - I m (\overset{̌}{f})$ [in (b) and (g), the dotted and solid lines overlap as in both cases, the gain parameter varies with time in exactly the same manner for paths AO and BO], while the third [(d) and (i)] and fourth [(e) and (j)] columns depict the temporal dynamics of the ratio factor r for two different evolution periods [κT = 30 (dark colored curves in the third column), κT′ = 5 (light colored curves in the fourth column), and T′ = T/6, green and orange lines denote the real and imaginary elements of r, respectively]. The filled and empty triangles (circles) at $\bar{t} = {\bar{t}}^{'} = 1$ in the last two columns describe the real and imaginary parts, respectively, of the component ratio associated with eigenstate (−) [(+)]—it can be clearly seen now how switching either to eigenmode (−) or (+) occurs for paths AO and BO, correspondingly, for sufficiently slowly varying conditions (column 3), thus confirming the analysis of Secs. II B 1 and II B 2. Initial excitation was randomly set, but kept the same between panels (d) and (e) panels (i) and (j) for a given parametric trajectory (AO or BO), to definitively illustrate that state conversion is an artifact of the magnitude of the total evolution period, rather than some other cause. In all instances, dotted/solid and dark/light lines indicate results relevant to paths AO or BO and periods T or T′, respectively, while the coupling constant was fixed in all simulation scenarios. It should also be noted that ɛ = 0.5 might not exactly satisfy the condition ɛ ≪ 1 used in Secs. II B 1 and II B 2 (this assumption was made to simplify our asymptotic calculations and to more efficiently convey the mathematical mechanism behind the observed mode switching effect) but was selected as such to demonstrate that state conversion emerges independently of the scale of ɛ.

The state vector dynamics can be retrieved after substituting f = η_f,1t + η_f,0 in Eq. , which takes the form

\frac{d^{2} a}{d t^{2}} + (- η_{f, 1}^{2} t^{2} - 2 η_{f, 1} η_{f, 0} t - η_{f, 0}^{2} + κ^{2} - η_{f, 1}) a = 0

(19)

and accepts as solutions the parabolic cylinder functions.³⁰ A first-order analysis with respect to ɛ enables us to omit term

η_{f, 1}^{2} t^{2}

from the preceding equation, which implies that Eq. becomes

\frac{d^{2} a}{d t^{2}} + (- 2 η_{f, 1} η_{f, 0} t - η_{f, 0}^{2} + κ^{2} - η_{f, 1}) a = 0 .

(20)

After applying the coordinate transformation,

χ = \frac{2 η_{f, 1} η_{f, 0} t + η_{f, 0}^{2} - κ^{2} + η_{f, 1}}{{(2 η_{f, 1} η_{f, 0})}^{2 / 3}} ≃ \frac{f^{2} - κ^{2} + η_{f, 1}}{h^{2}},

(21)

where

h = {(2 η_{f, 1} η_{f, 0})}^{1 / 3}

⁠, the standard Airy/Stokes equation is attained,

\frac{d^{2} a}{d χ^{2}} - χ a = 0,

(22)

which is satisfied by the Airy functions of the first a₁ = Ai(χ) and second kind a₂ = Bi(χ) (the reader should keep in mind that

f^{2} ≃ 2 η_{f, 1} η_{f, 0} t + η_{f, 0}^{2}

for small values of ɛ). The elements of the transfer matrix Z(T, 0) can now be directly computed from Eq. and the chain differentiation rule [da_1,2/dt = da_1,2/dχ, since

d χ / d t = h = {(2 η_{f, 1} η_{f, 0})}^{1 / 3}

], while it will also be true from Eq. that

W (t) = \frac{d χ}{d t} [A i (χ) \frac{d B i (χ)}{d χ} - B i (χ) \frac{d A i (χ)}{d χ}] = \frac{h}{π}

(23)

after utilizing the value of the Wronskian for the Airy functions.³¹

Of particular interest is to investigate the behavior of matrix Z(T, 0) in the slowly varying limit. In this regard and since

η_{f, 1} \overset{T \to \infty}{\to} 0

as can be seen from Eq. , it comes that

| χ | \overset{T \to \infty}{\to} \infty

(|·| stands for the complex absolute value function), and as such, appropriate Airy asymptotic expansions can be employed.³¹ For the trajectories depicted in Fig. 1(a), the subsequent relations apply (f₀ = κ ∓ iɛκ, f_T = κ − ɛκ) as follows:

f_{0}^{2} - κ^{2} = \mp 2 i ε κ^{2} (1 \mp i ε / 2) ≃ | f_{0}^{2} - κ^{2} | e^{\mp i π / 2 \mp i ε / 2},

(24a)

f_{T}^{2} - κ^{2} = κ^{2} (- 2 ε + ε^{2}) ≃ | f_{T}^{2} - κ^{2} | e^{i π},

(24b)

h = {(2 η_{f, 1} η_{f, 0})}^{1 / 3} ≃ {(\sqrt{2} ε κ^{2} T^{- 1})}^{1 / 3} e^{\pm i π / 4 \mp i ε / 3},

(24c)

along with

χ_{0} ≃ \frac{| f_{0}^{2} - κ^{2} |}{{(\sqrt{2} ε κ^{2} T^{- 1})}^{2 / 3}} e^{\mp i π \pm i ε / 6},

(25a)

χ_{T} ≃ \frac{| f_{T}^{2} - κ^{2} |}{{(\sqrt{2} ε κ^{2} T^{- 1})}^{2 / 3}} e^{\pm i π / 2 \pm i 2 ε / 3},

(25b)

where

| f_{0}^{2} - κ^{2} | ≃ | f_{T}^{2} - κ^{2} | ≃ 2 ε κ^{2}

and

| χ_{0} | ≃ | χ_{T} | ≃ {(2 \sqrt{ε} κ T)}^{2 / 3}

⁠. Given the angular sector in which the phases of χ₀, χ_T belong (−π to −2π/3, −2π/3 to 0, 0 to 2π/3, or 2π/3 to π), the corresponding Airy asymptotic formulas can be attained from Ref. 31, and after straightforward calculations, it can be found for very large periods T that

Z (T, 0) \sim [\begin{matrix} \frac{(e^{F_{T} - F_{0}} \pm i e^{F_{T} + F_{0}} \pm i e^{- F_{T} - F_{0}}) \sqrt[4]{χ_{0}}}{2 \sqrt[4]{χ_{T}}} & \frac{(e^{F_{T} - F_{0}} \mp i e^{F_{T} + F_{0}} \pm i e^{- F_{T} - F_{0}})}{2 h \sqrt[4]{χ_{0}} \sqrt[4]{χ_{T}}} \\ \frac{(e^{F_{T} - F_{0}} \pm i e^{F_{T} + F_{0}} \mp i e^{- F_{T} - F_{0}}) h \sqrt[4]{χ_{0}} \sqrt[4]{χ_{T}}}{2} & \frac{(e^{F_{T} - F_{0}} \mp i e^{F_{T} + F_{0}} \mp i e^{- F_{T} - F_{0}}) \sqrt[4]{χ_{T}}}{2 \sqrt[4]{χ_{0}}} \end{matrix}],

(26)

with

F_{0} = 2 / 3 χ_{0}^{3 / 2}

and

F_{T} = 2 / 3 χ_{T}^{3 / 2}

⁠. (The symbol ∼ denotes the asymptotic equivalence between functions, and it has an analogous interpretation as the big O notation. For details pertaining to the derivation of the Airy asymptotics, the reader can refer to Appendix A 3.) It can be seen now that the real part of the first (and dominating) term in the Puiseux series expansions of F₀ and F_T will be zero and negative

(- 2 \sqrt{2 ε} κ T / 3)

⁠, respectively. Consequently, in the limit T → ∞, it is expected that among the exponential factors appearing in the aforementioned expression for Z(T, 0), only the ones containing −F_T in the exponent will eventually survive. Hence, the asymptotic form of Z(T, 0) becomes

Z (T, 0) \sim [\begin{matrix} \frac{\pm i e^{- F_{T} - F_{0}} \sqrt[4]{χ_{0}}}{2 \sqrt[4]{χ_{T}}} & \frac{\pm i e^{- F_{T} - F_{0}}}{2 h \sqrt[4]{χ_{0}} \sqrt[4]{χ_{T}}} \\ \frac{\mp i e^{- F_{T} - F_{0}} h \sqrt[4]{χ_{0}} \sqrt[4]{χ_{T}}}{2} & \frac{\mp i e^{- F_{T} - F_{0}} \sqrt[4]{χ_{T}}}{2 \sqrt[4]{χ_{0}}} \end{matrix}] .

(27)

The following relations will hold for the elements of the above matrix:

\begin{align} \frac{ζ^{21} (T, 0)}{ζ^{11} (T, 0)} & ≃ \frac{ζ^{22} (T, 0)}{ζ^{12} (T, 0)} \to - h {\sqrt[4]{χ_{T}}}^{2} \\ ≃ \sqrt{| f_{T}^{2} - κ^{2} |} e^{\mp i π / 2} = \mp i \sqrt{κ^{2} - f_{T}^{2}}, \end{align}

(28)

which based on Eq. suggests that quantity (da/dt)/a will approach the value

\mp i \sqrt{κ^{2} - f_{T}^{2}}

for large evolution times. This in turn implies that the limiting behavior of the state vector component ratio can be retrieved from Eq. as

r_{T} \to \mp \sqrt{1 - \frac{f_{T}^{2}}{κ^{2}}} + i \frac{f_{T}}{κ} .

(29)

The right-hand side in the latter formula actually represents the ratio of the eigenvector components at time T, as computed in Appendix A 1 for

P T

-symmetric-like systems. More specifically, for the parametric paths AO (top sign) and BO (bottom sign) shown in Fig. 1(a), conversion to modes (−) and (+) is anticipated, respectively (state vector

|u_{T}⟩ = {[a_{T} b_{T}]}^{T}

and eigenvector

|v_{T}^{\mp}⟩ = {[v_{T, x}^{\mp} v_{T, y}^{\mp}]}^{T}

will essentially describe the same quantum state, since

r_{T} = b_{T} / a_{T} = v_{T, y}^{\mp} / v_{T, x}^{\mp}

⁠), and this is confirmed from the simulation results of Fig. 1(d) [in Fig. 1(e), it is shown that for smaller values of T, mode switching is not reached]. Meanwhile, if the ending point O of paths AO and BO lies to the right of the EP at (κ, 0), then after repeating the same calculations, conversion to the same eigenstate

|v_{T}^{-}⟩

was found. The delineated effect on the resulting system dynamics of the location of point O with respect to the EP degeneracy will be further explained in Sec. IV, in terms of the analysis provided in Sec. III and Appendix B 3, so as the nature of the aforementioned switching mechanism becomes fully comprehensible.

2. Hyperbolic time dependence for complex damping factor f(t)

Here, we shall examine the asymptotic properties of the

P T

-symmetric-like arrangement under study, when factor f varies in a hyperbolic fashion with time, i.e.,

f = {(η_{f, 1} t + η_{f, 0})}^{- 1}

⁠. Such a situation is illustrated in Figs. 1(f)–1(h) [A (κ, ɛκ), B (κ, −ɛκ), O (κ − ɛκ, 0)], with parameters η_f,1 and η_f,0 being evaluated as

η_{f, 1} = \frac{1}{κ T} (\frac{1}{1 - ε} - \frac{1}{1 \mp i ε}) ≃ \frac{\sqrt{2} ε}{κ T} e^{\mp i π / 4},

(30a)

η_{f, 0} = \frac{1}{κ \mp i ε κ} ≃ \frac{1}{κ} e^{\pm i ε} .

(30b)

In the preceding expressions, we assumed small values of

ε [e^{i ε} ≃ 1 + i ε, 1 / (1 + ε) ≃ 1 - ε, ε \in R^{+}]

⁠, while the top and bottom signs correspond to trajectories AO and BO in Fig. 1(f). Subsequently, we shall turn our attention to Eq. instead of Eq. , as it attains the more familiar Bessel-type form when

f = {(η_{f, 1} t + η_{f, 0})}^{- 1}

⁠. Nevertheless, to acquire the exact Bessel solutions, the change of variables =l(η_f,1t + η_f,0) and substitution

\tilde{a} = {(η_{f, 1} t + η_{f, 0})}^{p} a = {(t / l)}^{p} a

will be performed with p, l being both complex constants, as inspired by the Frobenius method. This results into Eq. being rewritten as

\frac{d^{2} a}{d t^{2}} + 2 (\frac{1}{η_{f, 1}} + p) \frac{d a}{d t} + (\frac{κ^{2}}{l^{2} η_{f, 1}^{2}} + \frac{p (p - 1) η_{f, 1} + 2 p}{t^{2} η_{f, 1}}) a = 0,

(31)

and after matching the coefficients with the modified Bessel differential equation, it will be true that p = −1/η_f,1 + 1/2 and

l^{2} = - κ^{2} / η_{f, 1}^{2}

⁠. Without loss of generality, we can choose l = iκ/η_f,1, and then, the general solution to Eq. is found to be

a = c_{1}^{a} I_{1 / η_{f, 1} - 1 / 2} (\frac{i κ}{f η_{f, 1}}) + c_{2}^{a} K_{1 / η_{f, 1} - 1 / 2} (\frac{i κ}{f η_{f, 1}}),

(32)

where I_ν(·) and K_ν(·) stand for the modified Bessel functions of the first and second kind of order ν, respectively. In terms of the original state vector component a, it will apply that

\begin{gathered} a = e^{\int_{0}^{t} f d t^{'}} \tilde{a} = e^{- \frac{\ln (f) - \ln (f_{0})}{η_{f_{1}}}} f^{- p} a = f^{- \frac{1}{2}} f_{0}^{\frac{1}{η_{f, 1}}} \\ a = c_{1}^{a} f^{- \frac{1}{2}} I_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}) + c_{2}^{a} f^{- \frac{1}{2}} K_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}), \end{gathered}

(33)

where χ = iκ/f,

c_{1, 2}^{a} = f_{0}^{1 / η_{f, 1}} c_{1, 2}^{a}

⁠, ln(·) is the natural logarithm function, and f₀ represents the value of factor f at t = 0. The full state vector dynamics can be directly retrieved from Eq. , with the elements of the transfer matrix Z(T, 0) given by Eq. after replacing

a_{1} = f^{- \frac{1}{2}} I_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}), a_{2} = f^{- \frac{1}{2}} K_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}),

(34)

and

\frac{d a_{1}}{d t} = \frac{1}{2} η_{f, 1} f a_{1} + i κ f^{- \frac{1}{2}} \frac{d I_{1 / η_{f, 1} - 1 / 2} (χ / η_{f, 1})}{d (χ / η_{f, 1})},

(35a)

\frac{d a_{2}}{d t} = \frac{1}{2} η_{f, 1} f a_{2} + i κ f^{- \frac{1}{2}} \frac{d K_{1 / η_{f, 1} - 1 / 2} (χ / η_{f, 1})}{d (χ / η_{f, 1})} .

(35b)

From the latter identities and Eq. , it can also be deduced that

\begin{align} W (t) & = i κ f^{- 1} [I_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}) \frac{d K_{1 / η_{f, 1} - 1 / 2} (χ / η_{f, 1})}{d (χ / η_{f, 1})} \\ - K_{1 / η_{f, 1} - 1 / 2} (\frac{χ}{η_{f, 1}}) \frac{d I_{1 / η_{f, 1} - 1 / 2} (χ / η_{f, 1})}{d (χ / η_{f, 1})}] = - η_{f, 1}, \end{align}

(36)

where we have utilized the value of the Wronskian for the modified Bessel functions.³² It should be noted that for the purposes of the current analysis, we could have employed the ordinary Bessel functions J_ν(·) (first kind) and Y_ν(·) (second kind) instead of I_ν(·) and K_ν(·). However, as it will become apparent from the subsequent adiabatic calculations, both I_ν(·) and K_ν(·) lead to simpler asymptotic expressions that involve exponentials, as opposed to J_ν(·) and Y_ν(·), which are described via the Airy functions in the asymptotic regime.³³

Having examined the analytic traits of the

P T

-symmetric-like configuration for the paths and temporal dependencies illustrated in Figs. 1(f)–1(h), we shall investigate the limiting behavior of the derived relations as T → ∞. From Eq. , it can be seen that

η_{f, 1} \overset{T \to \infty}{\to} 0

⁠, which in turn implies that the order of the modified Bessel functions governing the state vector dynamics can be approximated as ν = 1/η_f,1 − 1/2 ≃ 1/η_f,1. The succeeding expansions can now be deployed for large |ν|,³³

\begin{align} I_{ν} (ν χ) \sim \frac{e^{ν γ}}{\sqrt{2 π ν} \sqrt[4]{1 + χ^{2}}}, K_{ν} (ν χ) \sim \frac{\sqrt{π} e^{- ν γ}}{\sqrt{2 ν} \sqrt[4]{1 + χ^{2}}}, \end{align}

(37a)

\begin{align} \frac{d I_{ν} (ν χ)}{d (ν χ)} \sim \frac{e^{ν γ} \sqrt[4]{1 + χ^{2}}}{\sqrt{2 π ν} χ}, \frac{d K_{ν} (ν χ)}{d (ν χ)} \sim \frac{\sqrt{π} e^{- ν γ} \sqrt[4]{1 + χ^{2}}}{\sqrt{2 ν} χ}, \end{align}

(37b)

with γ defined as

γ = \sqrt{1 + χ^{2}} + \ln (\frac{χ}{1 + \sqrt{1 + χ^{2}}}),

(38)

\frac{d I_{ν} (ν χ) / d (ν χ)}{I_{ν} (ν χ)} ≃ \frac{d K_{ν} (ν χ) / d (ν χ)}{K_{ν} (ν χ)} \to G = \frac{\sqrt{1 + χ^{2}}}{χ} .

(39)

Of interest is to notice that applying the previous equation in Eq. and bearing in mind Eq. , we obtain

\frac{d a_{1, 2}}{d t} \sim (\frac{1}{2} η_{f, 1} f + i κ G) a_{1, 2} ≃ i κ G a_{1, 2},

(40)

since in the adiabatic limit, η_f,1 → 0. As already mentioned, Eqs. – are relevant when both ν and νγ lie in the right half of the complex plane. In any other case, appropriate analytic continuation formulas must be employed.³²

Let us proceed in evaluating the asymptotic form of Z(T, 0). Along these lines, it will be true that (f₀ = κ ∓ iɛκ, f_T = κ − ɛκ)

\arg (ν γ_{T}) ≃ \arg (ν γ_{0}) ≃ π / 2 \pm π / 4,

(41)

while

ν ≃ 1 / η_{f, 1} ≃ κ T e^{\pm i π / 4} / (\sqrt{2} ε)

I_{ν} (ν χ) = e^{i π ν^{'}} I_{ν^{'}} (ν^{'} χ^{'}),

(42a)

K_{ν} (ν χ) = - i π I_{ν^{'}} (ν^{'} χ^{'}) + e^{- i π ν^{'}} K_{ν^{'}} (ν^{'} χ^{'}),

(42b)

while for path BO, we simply have I_ν(νχ) = I_ν′(ν′χ′) and K_ν(νχ) = K_ν′(ν′χ′). For large periods T, Eqs. – can now be exploited, after replacing all quantities with their primed versions. [Substitution ν′ = ν might seem trivial above, but in many cases, ν may acquire a negative real part (like, for example, if point O was located to the right of the EP) and then we must set ν′ = −ν. The goal is to eventually have both ν′γ′ and ν′ lie in the right half of the complex plane, so as Eqs. – can be of use.] A few steps of direct computations will consequently lead to

Z (T, 0) \sim f_{T}^{- \frac{1}{2}} f_{0}^{- \frac{1}{2}} [\begin{matrix} \pm i κ G_{0}^{'} \cosh (ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}) & \sinh (ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}) \\ κ^{2} G_{T}^{'} G_{0}^{'} \sinh (ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}) & \pm i κ G_{T}^{'} \cosh (ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}) \end{matrix}],

(43)

where cosh(·) and sinh(·) represent the hyperbolic cosine and sine functions, respectively. In order to study the ratios ζ²¹(T, 0)/ζ¹¹(T, 0) and ζ²²(T, 0)/ζ¹²(T, 0), it is of essence to determine the asymptotic behavior of

ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}

⁠. The following Puiseux series expansions hold:

γ_{0}^{'} = \mp i \frac{π}{2} + \frac{2 \pm 2 i}{3} ε^{3 / 2} + O (ε^{5 / 2}),

(44a)

γ_{T}^{'} = \mp i \frac{π}{2} + \frac{2 i \sqrt{2}}{3} ε^{3 / 2} + O (ε^{5 / 2}),

(44b)

and thus

ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'} = \frac{κ T}{3} \sqrt{8 \mp 4 \sqrt{2}} e^{i 3 π / 4 \pm i 3 π / 8} ε^{1 / 2} + O (ε^{3 / 2}),

(45)

where O(·) signifies the big O notation. Since the real part of

ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}

is negative/positive for the segment AO/BO (top/bottom sign), it will apply that

\tanh (ν^{'} γ_{T}^{'} - ν^{'} γ_{0}^{'}) \to \mp 1

in the limit T → ∞, with tanh(·) designating the hyperbolic tangent [tanh(·) = sinh(·)/cosh(·)]. As such, we obtain

\frac{ζ^{21} (T, 0)}{ζ^{11} (T, 0)} ≃ \frac{ζ^{22} (T, 0)}{ζ^{12} (T, 0)} \to i κ G_{T}^{'} = \mp i κ G_{T} = \mp i \sqrt{κ^{2} - f_{T}^{2}},

(46)

where we have used the relation

G_{T}^{'} = \mp G_{T}

[given Eq. , it comes that

G^{'} = \sqrt{1 + χ^{' 2}} / χ^{'} = \mp \sqrt{1 + χ^{2}} / χ = \mp G

since χ′ = ∓χ]. Equation is equivalent to Eq. of Sec. II B 1, and as a result, the ratio of the state vector components at the end of the parametric evolution will obey Eq. . In other words, similarly to the examined scenarios in Sec. II B 1, switching to modes (−) and (+) is predicted for trajectories AO and BO, respectively, and this is confirmed by the numerical findings of Fig. 1(i) [for comparison purposes, analogous plots are provided in Fig. 1(j) but for a reduced evolution period]. If now point O lies to the right of the EP degeneracy, then the state vector is found to convert to eigenmode (−), in accordance with our calculations when factor f was assumed to depend linearly on time. The results of Secs. II B 1 and II B 2 essentially guide us to the subsequent key conclusion: for sufficiently slow driving and acyclic conditions, switching to one of the instantaneous eigenstates is expected at least for the analytically solvable traceless and symmetric non-Hermitian arrangement under study, which in turn constitutes a distinct violation of the conventional quantum adiabatic theorem. Moreover, we have illustrated from a mathematical standpoint, how such a dynamic trait is a consequence of the Stokes phenomenon of asymptotics, which can be summarized by the fact that the asymptotic expansions of certain classes of functions (here, the Airy and Bessel families) can vary in different regions of the complex plane. In what follows, we shall theoretically and numerically demonstrate how the observed state conversion effect is a standard feature of generalized (cyclic or acyclic) non-Hermitian Hamiltonian configurations (except for very special circumstances, which will be discussed), while also elaborating via particular examples on (i) the pivotal role of the eigenspectrum evolution graphs in determining the quantum state response and (ii) the existence of two separate but neighboring regimes of adiabatic rate values, which lead to differing switching dynamics.

III. ADIABATIC EVOLUTION PROPERTIES OF GENERALIZED HAMILTONIAN SETTINGS

Had we considered a time-independent Hamiltonian model [

M (t) = M_{o}, m^{j k} (t) = m_{o}^{j k}

with j, k = 1, 2], it becomes apparent from Eq. that the state vector component a(t) would obtain the following exact form:

a (t) = c_{a}^{+} e^{t λ_{o}^{+}} + c_{a}^{-} e^{t λ_{o}^{-}},

(47)

where

c_{a}^{\pm}

denote the complex amplitudes and

λ_{o}^{\pm}

are the roots to the characteristic polynomial associated with the second-order differential equation depicted in Eq. (or equivalently the roots to the characteristic polynomial of matrix M_o). If the Hamiltonian arrangement exhibits now a time dependence, which permits a perturbative treatment [M(t) = M_o + V(t), |V(t)/M_o| ≪ 1], then the solution to Eq. can be written as

a (t) = c_{a}^{+} (t) e^{t λ_{o}^{+}} + c_{a}^{-} (t) e^{t λ_{o}^{-}} .

(48)

The previous equation is similar to Eq. , except from amplitudes

c_{a}^{\pm} (t)

being time-dependent. Such a modification is necessary, since functions

e^{t λ_{o}^{\pm}}

will only approximately satisfy Eq. . [Had we retrieved the exact solutions a₁(t) and a₂(t), then the general (and exact) solution would take the form

a (t) = c_{1}^{a} a_{1} (t) + c_{2}^{a} a_{2} (t)

with

c_{1, 2}^{a}

being constants—see Sec. II A.] Moving even a step further, we will assume that the system matrix M(t) induces non-perturbative but adiabatic changes in the quantum state. Under such slowly varying conditions, it is reasonable to presume that the state vector dynamics will be still dictated by the inherent/natural quantum oscillation frequencies

[λ_{H}^{\pm} (t) = i λ^{\pm} (t)]

⁠, whose time dependence must be taken now into account. In this regard, the component a(t) can be expressed as

a (t) = c_{a}^{+} (t) e^{\int_{0}^{t} λ^{+} (t^{'}) d t^{'}} + c_{a}^{-} (t) e^{\int_{0}^{t} λ^{-} (t^{'}) d t^{'}},

(49)

which also falls in line with standard quantum adiabatic theory treatments [if

λ^{\pm} (t^{'}) = λ_{o}^{\pm}

⁠, then

\int_{0}^{t} λ^{\pm} (t^{'}) d t^{'} = t λ_{o}^{\pm}

and thus Eqs. and can be attained]. It is interesting to note that substituting

a (t) \to a^{\pm} (t) = e^{\int_{0}^{t} λ^{\pm} (t^{'}) d t^{'}}

into Eq. leads to the subsequent relation,

\begin{align} \frac{d λ^{\pm} (t)}{d t} + {[λ^{\pm} (t)]}^{2} - \{t r [M (t)] + \frac{d m^{12} (t)}{d t} m^{12} (t)\} λ^{\pm} (t) \\ + \det [M (t)] + m^{11} (t) / \frac{d m^{12} (t)}{d t} m^{12} (t) - \frac{d m^{11} (t)}{d t} = 0, \end{align}

(50)

which for slowly varying system parameters [

d Q (t) / d t \overset{T \to \infty}{≃} 0

⁠, with Q(t) = m^jk(t) or λ^±(t) and T denoting the overall period of evolution] becomes

{[λ^{\pm} (t)]}^{2} - t r [M (t)] λ^{\pm} (t) + \det [M (t)] = 0 .

(51)

The latter polynomial expression reflects the characteristic equation associated with matrix M(t) and is of course satisfied by the corresponding eigenvalues λ^±(t) at each time instant. In other words,

e^{\int_{0}^{t} λ^{\pm} (t^{'}) d t^{'}}

act as approximate and independent solutions to Eq. , and consequently, a(t) can be written as a linear superposition of such exponential factors at least in the adiabatic limit [Eq. ]. This has been more rigorously demonstrated via the Wentzel–Kramers–Brillouin (WKB) method in Appendix A 2, where explicit formulas for coefficients

c_{a}^{\pm} (t)

in Eq. have been retrieved. Along these lines, the results presented in Appendix A 2 can be employed to predict the asymptotic properties of a large class of special functions, frequently arising in the field of mathematical physics. We shall promptly illustrate this point for the Airy and modified Bessel functions, which emerged when analytically studying the

P T

-symmetric-like model in Sec. II B {

λ^{\pm} (t) = \pm i \sqrt{κ^{2} - f^{2} (t)}

from Eq. , since tr[M(t)] = 0 and det[M(t)] = κ² − f²(t) from Eq. }. In the former case, we have f(t) = η_f,1t + η_f,0 (Sec. II B 1), and thus, the following will apply [

f^{2} (t) ≃ 2 η_{f, 1} η_{f, 0} t + η_{f, 0}^{2}

for small ɛ]:

\begin{align} e^{\int_{0}^{t} λ^{\pm} (t^{'}) d t^{'}} & ≃ e^{\pm \int_{0}^{t} \sqrt{2 η_{f, 1} η_{f, 0} t^{'} + η_{f, 0}^{2} - κ^{2}} d t^{'}} \\ \propto e^{\pm \frac{2}{3} \frac{{(2 η_{f, 1} η_{f, 0} t + η_{f, 0}^{2} - κ^{2})}^{3 / 2}}{2 η_{f, 1} η_{f, 0}}} \\ ≃ e^{\pm \frac{2}{3} \frac{{[f^{2} (t) - κ^{2}]}^{3 / 2}}{h^{3}}} = e^{\pm \frac{2}{3} {[\frac{f^{2} (t) - κ^{2}}{h^{2}}]}^{3 / 2}} = e^{\pm \frac{2}{3} {[χ (t)]}^{3 / 2}}, \end{align}

(52)

with

h = {(2 η_{f, 1} η_{f, 0})}^{1 / 3}

⁠, χ(t) ≃ [f²(t) − κ²]/

h^{2}

for very large periods T [see Eq. , where

η_{f, 1} \overset{T \to \infty}{\to} 0

given Eq. ], and the proportionality symbol ∝ signifying an additional multiplicative constant stemming from the lower limit of integration. {Note that in Eq. [and in Eq. ], we have not fully accounted for the effect of principal roots, which may result in the exponent having an opposite sign.} The exponential elements in the above equation do capture the dynamics of Ai[χ(t)] and Bi[χ(t)] for |χ(t)| → ∞,³¹ but for a thorough derivation of the underlying asymptotics (including higher-order WKB terms and approximation errors), the reader should refer to Appendix A 3. An analogous computation can also be performed in the context of the modified Bessel functions. Given then that

f (t) = {(η_{f, 1} t + η_{f, 0})}^{- 1}

(Sec. II B 2) and after setting χ(t) = iκ/f(t) = iκ(η_f,1t + η_f,0) [χ(t) of Sec. II B 2 should not be confused with χ(t) of Sec. II B 1], the following relation will hold:

\begin{align} e^{\int_{0}^{t} λ^{\pm} (t^{'}) d t^{'}} & = e^{\pm i \int_{χ (0)}^{χ (t)} \sqrt{κ^{2} + \frac{κ^{2}}{χ^{2} (t^{'})}} \frac{d χ (t^{'})}{i κ η_{f, 1}}} \\ = e^{\pm \frac{1}{η_{f, 1}} \int_{χ (0)}^{χ (t)} \frac{\sqrt{1 + χ^{2} (t^{'})}}{χ (t^{'})} d χ (t^{'})} \\ \propto e^{\pm \frac{1}{η_{f, 1}} [\sqrt{1 + χ^{2} (t)} + \ln (\frac{χ (t)}{1 + \sqrt{1 + χ^{2} (t)}})]} ≃ e^{\pm ν γ (t)} . \end{align}

(53)

The preceding factors indeed outline the asymptotic form of I_ν[νχ(t)] and K_ν[νχ(t)],³³ with ν = 1/η_f,1 − 1/2 ≃ 1/η_f,1 for large T [

η_{f, 1} \overset{T \to \infty}{\to} 0

from Eq. ] and γ(t) being defined in Eq. (see Appendix A 4 for higher-order WKB terms and relative errors). It should be highlighted that in both Eqs. and , the absence of any spectrum discontinuities was assumed to simplify our calculations. Such an anomalous behavior can occur for instance after the evaluation of the complex square root function near the negative real axis, which may naturally come up in two-level Hamiltonian systems. (Higher-order complex roots will arise when considering higher-dimensional Hilbert spaces, but in all cases, the use of principal roots implies the emergence of branch cuts along the negative real axis—the impact of these discontinuity lines is discussed in detail within Appendix A 5 via specific examples.) Of course, the previous is more of a mathematical artifact, since from a physical perspective, eigenvalues and eigenvectors should be continuous functions of time, which is additionally a prerequisite for smoothly varying mode occupancy levels. In this vein, an alternative mode notation

[\tilde{(\pm)}]

is introduced in Fig. 2, whereby eigenvalues λ^±(t) are rearranged at each time instant so as to avoid any abrupt “jumps” in the eigenspectra evolution. Henceforth, the

\tilde{(\pm)}

mode representation shall be exclusively employed instead of (±), as it will (i) ensure the continuity of the involved time-dependent functions (and consequently the well-definedness of the corresponding time derivatives and integrals) and (ii) enable us to draw clear and useful conclusions based on the spectrum plots, which are not obscured by any discontinuous mode transitions.

FIG. 2.

View large Download slide

Eigenspectra evolution graphs for the $P T$ -symmetric-like model of Sec. II B. (a) and (b) Schematic of a rhombic parametric trajectory $C^{R}$ [(a)], along with the respective temporal dependencies of the gain (g) and detuning (σ) coefficients assuming clockwise encirclement in parameter space [(b)]. The starting point of evolution is signified with S, the center of the loop is located at (g_C, σ_C), while its horizontal and vertical extents are denoted by the real and positive variables ɛ_g, ɛ_σ. (c)–(f) Eigenspectra plots within the complex domain [ $R e (\overset{̌}{λ})$ ⁠, $I m (\overset{̌}{λ})$ ] corresponding to (±) [(c) and (e)] and $\tilde{(\pm)}$ [(d) and (f)] mode formalisms (⁠ $\overset{̌}{λ} = λ / κ, \overset{̌}{g} = g / κ, \overset{̌}{σ} = σ / κ, {\overset{̌}{ε}}_{g, σ} = ε_{g, σ} / κ$ ⁠, in which λ represents the eigenvalues of system matrix M). The notation $\tilde{(\pm)}$ refers to an appropriate reordering of modes (±) (⁠ $λ^{\pm} = \pm i \sqrt{κ^{2} - f^{2}}$ for the traceless and symmetric Hamiltonian model of Sec. II B), so as to achieve a smooth spectra evolution [as shown in (d) and (f)— $\tilde{(+)}$ is assigned to the mode branch, which at least originally emerges as dominant, i.e., $R e (λ^{\tilde{+}}) > R e (λ^{\tilde{-}})$ at t = 0], and avoid any discontinuities [such abrupt “jumps” appear in (c) and (e) owing to the evaluation of the principal square root in the formula for λ^±—for the definition of principal roots, see the beginning of Sec. II B 1—and are depicted in our graphs with dotted lines]. It can be deduced now from (d) that at the end of the parameter cycle, modes $\tilde{(-)}$ and $\tilde{(+)}$ are the instantaneous $(D)$ and averagely $(D_{a v})$ dominant eigenstates, respectively, while in (f), both $D$ and $D_{a v}$ are identified with $\tilde{(+)}$ ⁠. The parameter scenarios considered here are (i) ${\overset{̌}{g}}_{C} = 1, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = 1, {\overset{̌}{ε}}_{σ} = 0.5$ (middle row results) and (ii) ${\overset{̌}{g}}_{C} = 2, {\overset{̌}{σ}}_{C} = 0.5, {\overset{̌}{ε}}_{g} = 0.8, {\overset{̌}{ε}}_{σ} = 0.8$ (bottom row results). In all cases, the starting point S along $C^{R}$ is illustrated with a black filled circle (S^± and $S^{\tilde{\pm}}$ represent the respective points in the eigenspectra complex domain), while the EP degeneracy is signified with a magenta asterisk.

So far, we have demonstrated the potential of Eq. (49) in describing the limiting behavior of differential equations of a form analogous to Eq. (3), in alignment also with the WKB method analyzed in Appendix A 2 (the formulas and analysis can be directly extended to higher-order differential equations, ascribed to higher-dimensional Hamiltonian configurations). Of interest is now to investigate the circumstances where one of the two terms appearing in the right-hand side of Eq. (49) always becomes prevalent, independently of the initial excitation conditions. This physically implies that the quantum response will be dictated by one of the two instantaneous eigenmodes [either $\tilde{(+)}$ or $\tilde{(-)}$ ], and therefore, state conversion will have occurred. Evidently, the aforementioned scenario cannot emerge in Hermitian structures, since $λ^{\tilde{\pm}} (t)$ will be purely imaginary and thus both terms in Eq. (49) will be purely oscillatory and of a similar magnitude (see Appendix B 1). This is also in agreement with conventional quantum adiabatic theory, according to which if an arrangement is excited at one of its instantaneous eigenstates, then it remains in such an eigenstate as long as any external influences act on it slowly enough (the latter equivalently suggests that the response will depend on the initial conditions).

Next, let us examine the non-Hermitian case, where

λ^{\tilde{\pm}} (t)

are in general complex. It goes without saying that if

R e [λ^{\tilde{\pm}} (t)] = c (t)

holds for all t {e.g., this can happen if Re[M(t)] = c(t)I, with I being the identity matrix and

c (t) \in R

}, then standard adiabatic theory still applies after a simple shift of the respective eigenenergies. In other words, an appropriate imbalance must exist between the real components of

λ^{\tilde{\pm}} (t)

for any state conversion effects to appear, which can also be tentatively deduced from Eq. . More specifically, if on average

λ^{\tilde{\pm}} (t)

develops a larger real part than

λ^{\tilde{\mp}} (t)

⁠, i.e.,

R e {\int_{0}^{t} [λ^{\tilde{\pm}} (t^{'}) - λ^{\tilde{\mp}} (t^{'})] d t^{'}} > 0

⁠, then mode

\tilde{(\pm)}

is expected to dominate in the expression of a(t) for large periods of evolution T or alternatively slow adiabatic rates β (=1/T). Indeed, one can infer that

\begin{array}{l} |\frac{e^{\int_{0}^{t} λ^{\tilde{\pm}} (t^{'}) d t^{'}}}{e^{\int_{0}^{t} λ^{\tilde{\mp}} (t^{'}) d t^{'}}}| & = e^{Re \{\int_{0}^{t} [λ^{\tilde{\pm}} (t^{'}) - λ^{\tilde{\mp}} (t^{'})] d t^{'}\}} \\ = e^{T R e \{\int_{0}^{\bar{t}} [λ^{\tilde{\pm}} ({\bar{t}}^{'}) - λ^{\tilde{\mp}} ({\bar{t}}^{'})] d {\bar{t}}^{'}\}} \overset{T \to \infty}{\to} \infty, \end{array}

where |e^z| = e^Re(z) for a complex number z and

\int_{0}^{t} λ^{\tilde{\pm}} (t^{'}) d t^{'} = T \int_{0}^{\bar{t}} λ^{\tilde{\pm}} ({\bar{t}}^{'}) d {\bar{t}}^{'}

⁠, bearing in mind that the system matrix M(t) and its eigenvalues are explicit functions of the scaled time coordinate

\bar{t} = t / T

(see the beginning of Sec. II A). {For the sake of brevity, we shall be mostly using from now on

{⟨ λ^{\tilde{\pm}} ⟩}_{t}

in the place of

\int_{0}^{t} λ^{\tilde{\pm}} (t^{'}) d t^{'} / t

to signify the average of

λ^{\tilde{\pm}}

over [0, t].} Yet, such a calculation does not take into account the full effect of terms

c_{a}^{\tilde{\pm}} (t) e^{\int_{0}^{t} λ^{\tilde{\pm}} (t^{'}) d t^{'}}

⁠. This can be amended for, by simply substituting Eq. in Eq. and subsequently retrieving the differential equations satisfied by

c_{a}^{\tilde{\pm}} (t)

⁠. In this respect, the population dynamics associated with the modes of a generalized two-level non-Hermitian setting have been attained in Appendix B 1, while it has been shown (both theoretically and numerically) in Appendix B 3 that when superadiabatic rates are considered, it is not the averagely

(D_{a v})

but the instantaneous

(D)

dominant eigenstate, which governs the state vector evolution. {The former/latter refers to the mode branch leading to a larger value for the real part of

{⟨ λ^{\tilde{\pm}} ⟩}_{t}

/

λ^{\tilde{\pm}} (t)

(see the eigenspectra plots in Fig. 2, along with the caption clarifications). The special cases

R e ({⟨ λ^{\tilde{+}} ⟩}_{t}) = R e ({⟨ λ^{\tilde{-}} ⟩}_{t})

⁠,

R e [λ^{\tilde{+}} (t)] = R e [λ^{\tilde{-}} (t)]

⁠, and

λ^{\tilde{+}} (t) = λ^{\tilde{-}} (t)

(EPs) have all been addressed in Appendix B 3, while a very illustrative demonstration is also depicted in Fig. 6 of Sec. IV as to how to treat such scenarios.} Intuitively, this comes as no surprise since after extended temporal intervals, the state describing a dissipative physical system will lose track of both the initial conditions and the path

C^{R}

traversed in

R

space (a characteristic example is the classical damped harmonic oscillator, which, over long periods, eventually ceases motion and thus “forgets” its initial displacement/velocity and its intermediate trajectory). The previous is a consequence of the loss of energy/information within non-conservative environments, which is inextricably connected with the violation of Liouville’s theorem in classical (phase-space volume is not conserved over time) and quantum configurations (the total probability of finding the system in any state does not remain constant, commonly known as non-unitary evolution).

While in Appendix B 3, a comprehensive analytical treatment is provided pertaining to the quasistatic behavior of non-Hermitian Hamiltonian structures in terms of the mode population equations (with a parallel emphasis on the effect of the non-adiabatic mode couplings), such a study required a number of intricate mathematical manipulations and computations that would burden the reader with unnecessary details (and thus shift the attention from the key aspects of the paper) if included in the main text. For the remainder of this section, we shall demonstrate an alternative and more elegant approach, by examining the evolution traits of the ratio of the state vector components. The analysis will be conducted in the discrete time domain, which in turn will enable us to employ the rich properties of Möbius transformations (appropriately adapted for the case of time-dependent Hamiltonians) in order to predict the dynamics of generalized two-mode systems under slowly varying conditions (the results can be straightforwardly extended to continuous time arrangements, in the limit of infinitesimally small time steps). In that regard, we shall direct our attention to Eq. , which can be discretized as follows (see also Appendix D):

|u_{n + 1}⟩ = [\begin{matrix} 1 + m_{n}^{11} Δ & m_{n}^{12} Δ \\ m_{n}^{21} Δ & 1 + m_{n}^{22} Δ \end{matrix}] |u_{n}⟩ = Ξ_{n} |u_{n}⟩ .

(54)

In the preceding expression,

n \in N \land n \in [0, N]

(T = NΔ signifies the overall evolution period),

|u_{n}⟩ = {[a_{n} b_{n}]}^{T}

⁠, Δ corresponds to the small and positive discrete time step, q_n = q(t = t_n = nΔ) stands for a discrete and time-dependent quantity, and Ξ_n = I + M_nΔ depicts the propagator matrix with elements

ξ_{n}^{j k} = δ^{j k} + m_{n}^{j k} Δ

(j, k = 1, 2, δ^jk: Kronecker delta symbol). The ratio factor r_n = b_n/a_n will now satisfy the following relation:

r_{n + 1} = G_{n} (r_{n}) = \frac{ξ_{n}^{21} + ξ_{n}^{22} r_{n}}{ξ_{n}^{11} + ξ_{n}^{12} r_{n}},

(55)

which represents a time-dependent Möbius transformation [in its continuous time manifestation, the above formula delineates a nonlinear differential equation of the Riccati-type, as specified in Eq. of Sec. II A—see Appendix D for more details]. Equation in essence maps r_n to r_n+1 based on the elements of Ξ_n—we shall be referring to such a mapping as the

G

-transform (its mathematical properties are discussed in Appendix E). A special class of points, known as “fixed points,” exist for this family of transforms, which satisfy the identity

G_{n} (r_{n}^{\tilde{\pm}}) = r_{n}^{\tilde{\pm}}

(the term “fixed points” is more of a misnomer here, since their position in the complex plane varies with n owing to the temporal dependence of

ξ_{n}^{j k}, m_{n}^{j k}

⁠, but we shall adopt the traditional terminology). It will be true then that

{(r_{n}^{\tilde{\pm}})}^{2} - \frac{m_{n}^{22} - m_{n}^{11}}{m_{n}^{12}} r_{n}^{\tilde{\pm}} - \frac{m_{n}^{21}}{m_{n}^{12}} = 0 .

(56)

Let us assume the following solutions to the aforementioned quadratic formula:

r_{n}^{\tilde{\pm}} = \frac{1}{m_{n}^{12}} λ_{n}^{\tilde{\pm}} - \frac{m_{n}^{11}}{m_{n}^{12}},

(57)

where

λ_{n}^{\tilde{\pm}} = λ^{\tilde{\pm}} (t = t_{n} = n Δ)

are the characteristic roots to Eq. . Such a postulation was not arbitrary, but instead was attained after employing Eq. of Sec. II A in its discrete form [an analogous expression can also be derived from Eq. ] and applying it for eigenstates

\tilde{(\pm)}

(r_{n} \to r_{n}^{\tilde{\pm}}, a_{n} \to a_{n}^{\tilde{\pm}})

⁠, as explicitly shown in Appendix D [in this sense,

r_{n}^{\tilde{\pm}}

denote the instantaneous eigenvector component ratios, which can be further corroborated via the eigenequation

(λ_{n}^{\tilde{\pm}} I - M_{n}) |v_{n}^{\tilde{\pm}}⟩ = 0

]. To validate our guess, we will resort to Vieta’s identities associated with Eq. , i.e.,

λ_{n}^{\tilde{+}} + λ_{n}^{\tilde{-}} = t r (M_{n}), λ_{n}^{\tilde{+}} \cdot λ_{n}^{\tilde{-}} = \det (M_{n})

⁠. Given then Eq. , it comes that

r_{n}^{\tilde{+}} + r_{n}^{\tilde{-}} = (m_{n}^{22} - m_{n}^{11}) / m_{n}^{12}

and

r_{n}^{\tilde{+}} \cdot r_{n}^{\tilde{-}} = - m_{n}^{21} / m_{n}^{12}

⁠, which correspond to Vieta’s formulas of Eq. . Consequently, Eq. accurately describes the solutions to Eq. .

It would be of particular interest at this moment, to physically interpret our choice to correlate $r_{n}^{\tilde{\pm}}$ (i.e., the “fixed points” of $G_{n}$ ⁠) with eigenvectors $|v_{n}^{\tilde{\pm}}⟩$ ⁠. In this vein, we can consider the relation $G_{n} (r_{n}^{\tilde{\pm}}) = r_{n}^{\tilde{\pm}}$ in conjunction with the conventional adiabatic theorem underpinning Hermitian quantum processes: if the instantaneous eigenstates are originally excited $(r_{0} = r_{0}^{\tilde{\pm}})$ ⁠, then these will continue to dictate the state vector dynamics in the slowly varying limit $[G_{n + 1} (G_{n} (r_{n}^{\tilde{\pm}})) = G_{n + 1} (r_{n}^{\tilde{\pm}}) ≃ G_{n + 1} (r_{n + 1}^{\tilde{\pm}}) = r_{n + 1}^{\tilde{\pm}}]$ ⁠. In the context of non-Hermitian settings, the latter statement will not uphold, since $r_{n}^{\tilde{\pm}}$ are not any more indifferent (i.e., they become attractive/repulsive—see the subsequent analysis), and therefore, small variations in the arguments of the $G$ -transform can cumulatively lead to large deviations in the output.

Following the standard theory of Möbius transformations, we evaluate the derivative of

G_{n} (r)

at

r_{n}^{\tilde{\pm}}

⁠,

{\frac{d G_{n}}{d r}|}_{r_{n}^{\tilde{\pm}}} = \frac{\det (Ξ_{n})}{{(ξ_{n}^{11} + ξ_{n}^{12} r_{n}^{\tilde{\pm}})}^{2}} .

(58)

The numerator in the above formula can be equivalently expressed as

\begin{align} \det (Ξ_{n}) & = ξ_{n}^{11} ξ_{n}^{22} - ξ_{n}^{12} ξ_{n}^{21} \\ = (1 + m_{n}^{11} Δ) (1 + m_{n}^{22} Δ) - m_{n}^{12} m_{n}^{21} Δ^{2} \\ = 1 + (m_{n}^{11} + m_{n}^{22}) Δ + (m_{n}^{11} m_{n}^{22} - m_{n}^{12} m_{n}^{21}) Δ^{2} \\ = [Δ^{- 2} + t r (M_{n}) Δ^{- 1} + \det (M_{n})] Δ^{2} = F (- Δ^{- 1}) Δ^{2}, \end{align}

(59)

where

F (z) = z^{2} - t r (M_{n}) z + \det (M_{n}) = (z - λ_{n}^{\tilde{+}}) (z - λ_{n}^{\tilde{-}})

designates the characteristic polynomial of matrix M_n with roots

λ_{n}^{\tilde{\pm}}

⁠. Moreover, the denominator in Eq. becomes equal to [see also Eq. ]

ξ_{n}^{11} + ξ_{n}^{12} r_{n}^{\tilde{\pm}} = 1 + m_{n}^{11} Δ + m_{n}^{12} \frac{λ_{n}^{\tilde{\pm}} - m_{n}^{11}}{m_{n}^{12}} Δ = 1 + λ_{n}^{\tilde{\pm}} Δ .

(60)

Based now on Eqs. and , and the form of the characteristic polynomial

F (z)

⁠, Eq. can be rewritten as

\begin{align} {\frac{d G_{n}}{d r}|}_{r_{n}^{\tilde{\pm}}} & = \frac{F (- Δ^{- 1}) Δ^{2}}{{(1 + λ_{n}^{\tilde{\pm}} Δ)}^{2}} = \frac{(- Δ^{- 1} - λ_{n}^{\tilde{+}}) (- Δ^{- 1} - λ_{n}^{\tilde{-}}) Δ^{2}}{{(1 + λ_{n}^{\tilde{\pm}} Δ)}^{2}} \\ = \frac{(1 + λ_{n}^{\tilde{+}} Δ) (1 + λ_{n}^{\tilde{-}} Δ)}{{(1 + λ_{n}^{\tilde{\pm}} Δ)}^{2}} ≃ e^{(λ_{n}^{\tilde{\mp}} - λ_{n}^{\tilde{\pm}}) Δ}, \end{align}

(61)

where we have taken into account that

e^{λ_{n}^{\tilde{\pm}} Δ} ≃ 1 + λ_{n}^{\tilde{\pm}} Δ

for small Δ. Comparing the absolute value of the expression in Eq. to unity will allow us to identify the type of

G

-transform along with its effect on the dynamics of r_n. In this respect, the succeeding scenarios can be discerned as follows:

Under Hermitian conditions $[\forall n : R e (λ_{n}^{\tilde{\pm}}) = 0]$ ⁠, $r_{n}^{\tilde{\pm}}$ are neither attractive nor repulsive but indifferent and the $G$ -transform can be classified as elliptic.
Under non-Hermitian driving, if $R e (λ_{n}^{\tilde{+}}) = R e (λ_{n}^{\tilde{-}})$ for all possible values of n, then the conclusions in (i) are also relevant here.
Under non-Hermitian driving [which in general implies $R e (λ_{n}^{\tilde{-}}) \neq R e (λ_{n}^{\tilde{+}})$ ], if $R e (λ_{n}^{\tilde{-}}) > R e (λ_{n}^{\tilde{+}})$ ⁠, then $r_{n}^{\tilde{+}}$ $(r_{n}^{\tilde{-}})$ is repulsive (attractive) and the transform is said to be loxodromic [if $R e (λ_{n}^{-}) < R e (λ_{n}^{+})$ ⁠, the roles of $r_{n}^{\tilde{+}}, r_{n}^{\tilde{-}}$ are reversed—in Appendix E, it is analytically shown how the attractive/repulsive behavior of $r_{n}^{\tilde{+}}, r_{n}^{\tilde{-}}$ emerges in loxodromic transformations].

From the preceding observations, we can infer that any initially excited quantum state will convert to the instantaneous dominant eigenmode $D$ [ $\tilde{(\pm)}$ ⁠, if $R e (λ_{n}^{\tilde{\pm}}) > R e (λ_{n}^{\tilde{\mp}})$ ] during adiabatic evolution [or mathematically stated, the iterative action of the $G$ -transform will result in $r_{n + 1} = (G_{n} ◦ G_{n - 1} ◦ \dots ◦ G_{0}) (r_{0}) \overset{n, N \to \infty}{\to} r_{n + 1}^{\tilde{\pm}}$ ⁠, with symbol ◦ indicating function composition and r₀ given by initial conditions], as long as non-Hermiticity and an imbalance between the real parts of $λ_{n}^{\tilde{\pm}}$ are both present. This in turn confirms the detailed investigation of Appendix B 3, at least when superadiabatic rates β are considered. Moreover, the impact of the various types of $G$ -transform on the evolution traits of r_n has been numerically demonstrated in Fig. 3 and is in agreement with our theoretical predictions. Finally, the reader can refer to the supplementary material for a generalization of our approach to three-level Hamiltonian systems, whereby the properties of the two-dimensional $G$ -transform have been employed.

FIG. 3.

View large Download slide

Discrete evolution dynamics of the ratio factor r_n according to the time-dependent Möbius (or $G$ -) transform in the adiabatic limit. The numerical results demonstrated here [(b)–(e)] correspond to the parametric trajectories shown in panel (a), with both gain and detuning factors depending linearly on time as analyzed in Sec. II B 1 and illustrated in Figs. 1(a)–1(c) [ ${\overset{̌}{f}}_{0} = 1 - i ε s$ and ${\overset{̌}{f}}_{N} = 1 - ε$ for the brown (s = 1), purple (s = 0.2), and gray (s = 0) paths; ${\overset{̌}{f}}_{0} = - i ε$ and ${\overset{̌}{f}}_{N} = 0$ for the light blue path; N = 500, κΔ = 0.1, ɛ = 0.5 in all cases]. The system matrix M_n takes the form shown in Eq. (14) as is [(b)–(d)] and shifted also by q_nI [(e)], i.e., $M_{n}^{'} = M_{n} + q_{n} I$ (I: 2 × 2 identity matrix; q_n: arbitrary real shift). The red and blue lines in (b)–(e) represent the component ratios of the instantaneous eigenvectors $|v_{n}^{\tilde{+}}⟩$ $(r_{n}^{\tilde{+}})$ and $|v_{n}^{\tilde{-}}⟩$ $(r_{n}^{\tilde{-}})$ ⁠, respectively, whereas the brown, purple, gray, and light blue scatter diagrams describe the dynamics of r_n associated with the same colored parametric trajectories in (a) [main panels and insets depict the evolution within the complex domain (arrows indicate the temporal direction) and with time (the absolute value of r_n is plotted against ${\bar{t}}_{n} = t_{n} / T = n / N$ ⁠), respectively—in (d) and (e), magnified graphs are additionally provided for a better illustration of the $r_{n}^{\tilde{\pm}}$ curves]. It is clear that (b) and (c) delineate a loxodromic $G$ -transform (the $r_{n}^{\tilde{+}}$ curve acts as an attractor for the r_n evolution line), where conversion to eigenstate $\tilde{(+)}$ takes place. This confirms the findings of Sec. II B 1 and Figs. 1(a)–1(e) [the correspondence between notations $\tilde{(\pm)}$ and (±) is shown in the spectrum plots of Fig. 6(a)], while it can also be deduced from (b) and (c) (and their insets) that the greater the absolute separation $| R e ({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) |$ [for the $P T$ -symmetric-like system under study, the previous quantity becomes equal to $2 | R e (i \sqrt{1 - {\overset{̌}{f}}_{n}^{2}}) | = 2 | I m (\sqrt{1 - {\overset{̌}{f}}_{n}^{2}}) |$ and attains overall a larger value for the brown as opposed to the purple path in (a), since the former is further distant apart from the EP located at $({\overset{̌}{g}}_{n}, {\overset{̌}{σ}}_{n}) = (1, 0)$ (magenta asterisk in the top left panel) compared to the latter], the faster the state vector converts to the dominant eigenmode [here $\tilde{(+)}$ ]. To the same conclusion we can arrive given Eq. (61), according to which the magnitude of factor $e^{({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) κ Δ}$ (bear in mind the relation $λ_{n}^{\tilde{\pm}} = {\overset{̌}{λ}}_{n}^{\tilde{\pm}} κ$ ⁠) and its comparison with unity determines the resulting degree of “loxodromicity.” Finally, panels (d) [ ${\overset{̌}{g}}_{n} \neq 0, {\overset{̌}{σ}}_{n} = 0, R e ({\overset{̌}{λ}}_{n}^{\tilde{\pm}}) = 0$ ] and (e) [ ${\overset{̌}{g}}_{n} = 0, {\overset{̌}{σ}}_{n} \neq 0, R e ({\overset{̌}{λ}}_{n}^{\tilde{\pm}}) = {\overset{̌}{q}}_{n} = q_{n} / κ$ —see also the temporal evolution of ${\overset{̌}{λ}}_{n}^{\tilde{\pm}}$ in (f), where the green and orange lines signify the real and imaginary parts, respectively, for modes $\tilde{(+)} / \tilde{(-)}$ (dark/light colored curves), while the inset indicates the ${\overset{̌}{λ}}_{n}^{\tilde{\pm}}$ plots when no shift is applied to M_n] belong to the class of elliptic transforms and exhibit zero degree of “loxodromicity” (⁠ $| e^{({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) κ Δ} | = 1$ ⁠). This is reflected by the fact that r_n no longer becomes attracted (repulsed) by the $r_{n}^{\tilde{\pm}}$ lines but instead tends to move in circles around them.

FIG. 3.

View large Download slide

Discrete evolution dynamics of the ratio factor r_n according to the time-dependent Möbius (or $G$ -) transform in the adiabatic limit. The numerical results demonstrated here [(b)–(e)] correspond to the parametric trajectories shown in panel (a), with both gain and detuning factors depending linearly on time as analyzed in Sec. II B 1 and illustrated in Figs. 1(a)–1(c) [ ${\overset{̌}{f}}_{0} = 1 - i ε s$ and ${\overset{̌}{f}}_{N} = 1 - ε$ for the brown (s = 1), purple (s = 0.2), and gray (s = 0) paths; ${\overset{̌}{f}}_{0} = - i ε$ and ${\overset{̌}{f}}_{N} = 0$ for the light blue path; N = 500, κΔ = 0.1, ɛ = 0.5 in all cases]. The system matrix M_n takes the form shown in Eq. (14) as is [(b)–(d)] and shifted also by q_nI [(e)], i.e., $M_{n}^{'} = M_{n} + q_{n} I$ (I: 2 × 2 identity matrix; q_n: arbitrary real shift). The red and blue lines in (b)–(e) represent the component ratios of the instantaneous eigenvectors $|v_{n}^{\tilde{+}}⟩$ $(r_{n}^{\tilde{+}})$ and $|v_{n}^{\tilde{-}}⟩$ $(r_{n}^{\tilde{-}})$ ⁠, respectively, whereas the brown, purple, gray, and light blue scatter diagrams describe the dynamics of r_n associated with the same colored parametric trajectories in (a) [main panels and insets depict the evolution within the complex domain (arrows indicate the temporal direction) and with time (the absolute value of r_n is plotted against ${\bar{t}}_{n} = t_{n} / T = n / N$ ⁠), respectively—in (d) and (e), magnified graphs are additionally provided for a better illustration of the $r_{n}^{\tilde{\pm}}$ curves]. It is clear that (b) and (c) delineate a loxodromic $G$ -transform (the $r_{n}^{\tilde{+}}$ curve acts as an attractor for the r_n evolution line), where conversion to eigenstate $\tilde{(+)}$ takes place. This confirms the findings of Sec. II B 1 and Figs. 1(a)–1(e) [the correspondence between notations $\tilde{(\pm)}$ and (±) is shown in the spectrum plots of Fig. 6(a)], while it can also be deduced from (b) and (c) (and their insets) that the greater the absolute separation $| R e ({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) |$ [for the $P T$ -symmetric-like system under study, the previous quantity becomes equal to $2 | R e (i \sqrt{1 - {\overset{̌}{f}}_{n}^{2}}) | = 2 | I m (\sqrt{1 - {\overset{̌}{f}}_{n}^{2}}) |$ and attains overall a larger value for the brown as opposed to the purple path in (a), since the former is further distant apart from the EP located at $({\overset{̌}{g}}_{n}, {\overset{̌}{σ}}_{n}) = (1, 0)$ (magenta asterisk in the top left panel) compared to the latter], the faster the state vector converts to the dominant eigenmode [here $\tilde{(+)}$ ]. To the same conclusion we can arrive given Eq. (61), according to which the magnitude of factor $e^{({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) κ Δ}$ (bear in mind the relation $λ_{n}^{\tilde{\pm}} = {\overset{̌}{λ}}_{n}^{\tilde{\pm}} κ$ ⁠) and its comparison with unity determines the resulting degree of “loxodromicity.” Finally, panels (d) [ ${\overset{̌}{g}}_{n} \neq 0, {\overset{̌}{σ}}_{n} = 0, R e ({\overset{̌}{λ}}_{n}^{\tilde{\pm}}) = 0$ ] and (e) [ ${\overset{̌}{g}}_{n} = 0, {\overset{̌}{σ}}_{n} \neq 0, R e ({\overset{̌}{λ}}_{n}^{\tilde{\pm}}) = {\overset{̌}{q}}_{n} = q_{n} / κ$ —see also the temporal evolution of ${\overset{̌}{λ}}_{n}^{\tilde{\pm}}$ in (f), where the green and orange lines signify the real and imaginary parts, respectively, for modes $\tilde{(+)} / \tilde{(-)}$ (dark/light colored curves), while the inset indicates the ${\overset{̌}{λ}}_{n}^{\tilde{\pm}}$ plots when no shift is applied to M_n] belong to the class of elliptic transforms and exhibit zero degree of “loxodromicity” (⁠ $| e^{({\overset{̌}{λ}}_{n}^{\tilde{-}} - {\overset{̌}{λ}}_{n}^{\tilde{+}}) κ Δ} | = 1$ ⁠). This is reflected by the fact that r_n no longer becomes attracted (repulsed) by the $r_{n}^{\tilde{\pm}}$ lines but instead tends to move in circles around them.

IV. SYMMETRIC AND ASYMMETRIC STATE CONVERSION UNDER THE ACTION OF CYCLIC AND SLOWLY VARYING NON-HERMITIAN HAMILTONIANS

Having provided an extensive qualitative and quantitative insight as to the quasistatic response of two-mode non-Hermitian quantum configurations traversing trajectories of any shape in $R$ space, we are now in the position to apply our findings in the case of closed paths $C^{R}$ ⁠. Before doing so, it is necessary to first specify the meaning of symmetric (SSC) and asymmetric (ASC) state conversion in the context of cyclic Hamiltonian arrangements: in the former/latter instance, clockwise (CW) and counterclockwise (CCW) encirclements within $R$ space lead to switching to the same/different instantaneous eigenstates after a full cycle, independently of how the system was initially prepared. Along these lines, and leveraging on the analysis of Sec. III and Appendix B 3, it can be deduced that for superadiabatic parameter changes, which entail conversion to mode $D$ ⁠, SSC will always be observed if $R e [λ^{\tilde{+}} (T)] \neq R e [λ^{\tilde{-}} (T)]$ (T: cyclic evolution period)—this becomes self-evident, since the eigenstate exhibiting a larger value for Re[λ(T)] (and therefore recognized as the $D$ mode) is uniquely determined and will not depend on the encircling direction. The only situation where ASC might possibly arise is when $R e [λ^{\tilde{+}} (T)] = R e [λ^{\tilde{-}} (T)]$ ⁠, but then $D$ cannot be strictly defined at t = T. This can be resolved by simply inspecting, just before time T, which of the $\tilde{(\pm)}$ eigenmode branches is characterized by instantaneous eigenvalues with a larger real part—such a branch will eventually dominate for sufficiently large T. The aforementioned rationale mathematically implies that the dynamical behavior of the quantum system is examined at the limit t → T, which in turn based on the theory of limits and the continuity of the underlying functions yields an accurate description of what happens at t = T. {For a graphical representation of the latter, we direct the reader to Figs. 6(a) and 6(b), where the dynamics along SE [ $R e (λ^{\tilde{\pm}}) = 0$ at E] are explored in terms of path SK, as point K approaches E.} A completely analogous reasoning can also be employed in the presence of eigenspectrum degeneracies, whether they occur at the end or anywhere along $C^{R}$ [see Figs. 6(c)–6(h)].

To substantiate our discussion thus far, separate numerical calculations have been performed, and the results are shown in Figs. 4–6 for the non-Hermitian model of Sec. II B, assuming rhombic parametric trajectories $C^{R}$ ⁠. In this vein, ASC was noticed in the first and third rows of Fig. 4 [see Figs. 4(b) and 4(f)], where the ending point of $C^{R}$ is located to the left of the EP degeneracy at (κ, 0), thereby leading to $R e [λ^{\tilde{\pm}} (T)] = 0$ ⁠. {The $P T$ -symmetric-like system under study exhibits eigenvalues $λ^{\pm} (t) = \pm i \sqrt{κ^{2} - f^{2} (t)}$ ⁠. Therefore, the relation $R e [λ^{\tilde{+}} (T)] = R e [λ^{\tilde{-}} (T)]$ can be satisfied iff $R e [λ^{\tilde{\pm}} (T)] = 0 \Leftrightarrow | g (T) | ⩽ κ, σ (T) = 0$ ⁠.} Meanwhile, in the second row of Fig. 4 and for all the simulation scenarios in Fig. 5, SSC was attained as then $R e [λ^{\tilde{+}} (T)] \neq R e [λ^{\tilde{-}} (T)]$ ⁠. (Here, we refer to the main panels of the last column of Fig. 5, where conversion to the $D$ mode takes place.) It should be highlighted that in all the aforementioned examples, switching to mode $D$ is observed for both CW and CCW parametric encirclements. This aligns with the analytical investigations of Sec. II B, where the quantum state was found to be dictated by the $D$ mode branch in the extreme superadiabatic limit T → ∞ [the findings of Sec. II B have been numerically verified through Figs. 1, 3(a)–3(c), 6(a), and 6(b), bearing also in mind that $D = \tilde{(+)} = (-)$ in these instances, as detailed in the caption of Fig. 6]. Such a dynamical behavior persists, even when EPs are crossed by $C^{R}$ —this is illustrated in Figs. 6(c)–6(h) [see the inset of Fig. 6(f)], with the state vector dynamics along $C^{R}$ being examined in terms of a family of auxiliary trajectories $C^{R}$ ⁠, which do not pass through any spectrum degeneracies (see the caption of Fig. 6, for a more comprehensive description).

FIG. 4.

View large Download slide

EP encirclement as not a necessary or sufficient requirement for ASC to occur. The top insets in the first column depict the paths $C^{R}$ within $R$ space [the rhombic geometry is utilized for our numerical calculations (see Fig. 2), S: starting point along $C^{R}$ ⁠, magenta asterisk: EP], while the main panels and bottom insets (which share the same axis margins) illustrate the $\tilde{(+)}$ (red) and $\tilde{(-)}$ (blue) eigenmode branches [ $S^{\tilde{\pm}}$ denote the starting points along the $\tilde{(\pm)}$ branches and the averaged quantity $⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩$ over period T is indicated with a partially filled red/blue circle outlined with a gray color] under CW and CCW conditions (CW/CCW related results are signified with solid/dotted lines in all graphs). In the second column, the green/orange curves and filled/empty circles or triangles describe the real/imaginary parts of the ratio factor r and of the eigenvector component ratios $r^{\tilde{\pm}}$ at the end of the cycle, respectively (circles: $r_{CW}^{\tilde{+}}$ ⁠; triangles: $r_{CCW}^{\tilde{-}}$ ⁠). (a) and (b) Parameters used: ${\overset{̌}{g}}_{C} = 1, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = 0.1, {\overset{̌}{ε}}_{σ} = 0.3, κ T = 100$ ⁠. The dominant $\tilde{(+)}$ mode branch $[D = D_{a v} = \tilde{(+)}]$ ends at different points for CW and CCW evolution, thus leading to asymmetric state conversion as confirmed by (b). (c) and (d) Parameters are identical to the first row, with the only difference being the starting point along $C^{R}$ ⁠. Symmetric mode switching is noticed from (d), despite $C^{R}$ enclosing the EP degeneracy, and this is attributed to the prevalent $\tilde{(-)}$ branch $[D = \tilde{(-)}, D_{a v} = \tilde{(+)}]$ terminating at the same location independently of the encircling direction. (e) and (f) Parameters employed: ${\overset{̌}{g}}_{C} = 0.6, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = {\overset{̌}{ε}}_{σ} = 0.3, κ T = 100$ ⁠. The graph in (f) shows that ASC occurs even for a non-EP enclosing loop, which can also be inferred from the differing ending points of the dominant $\tilde{(-)}$ branch $[D = \tilde{(-)}, D_{a v} = \tilde{(+)}]$ under CW and CCW parametric encirclements. In all instances, conversion to the $D$ eigenstate takes place due to superadiabatic conditions (this is trivial for the first row, where $D = D_{a v}$ ⁠), while the assignment of $D, D_{a v}$ to modes $\tilde{(\pm)}$ comes after comparing $R e ({\overset{̌}{λ}}^{\tilde{\pm}})$ and $R e (⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩)$ ⁠, respectively, near or at $\bar{t} = t / T = 1$ ⁠.

FIG. 5.

View large Download slide

Reconfigurable operation between the ASC and SSC switching regimes. The first two columns illustrate the closed trajectories $C^{R}$ within $R$ space (insets) along with the evolution of the eigenvalues ${\overset{̌}{λ}}^{\tilde{\pm}}$ in the complex domain (main panels), under clockwise (first column) and counterclockwise (second column) conditions. The last two columns depict the dynamics of the ratio factor for moderate (third column—T₁) and large (fourth column—T₂) values of the cycle period T (state conversion occurs in both occasions, but to modes $D_{a v}$ and $D$ ⁠, respectively, as we shall see). The colors and types of lines and symbols appearing in the various graphs have already been described in Fig. 4. (a)–(d) Trajectory parameters: ${\overset{̌}{g}}_{C} = 1, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = {\overset{̌}{ε}}_{σ} = 0.2, \bar{τ} = 1 / 16, κ T_{1} = 100, κ T_{2} = 2500$ ⁠. [Variable $\bar{τ} = τ / T$ denotes the CW $(\bar{τ} > 0)$ or CCW $(\bar{τ} < 0)$ normalized time shift from the left corner of $C^{R}$ ⁠.] For CW encircling [(a)], it can be inferred that $D_{a v} = \tilde{(+)}$ and $D = \tilde{(-)}$ ⁠, while in the CCW scenario [(b)], $D = D_{a v} = \tilde{(-)}$ ⁠. Consequently, and according to the analysis of Sec. III and Appendix B 3, conversion to the $D_{a v}$ mode is originally expected. Yet, as T is increased, there will be a critical transition point (T_cr) where $D_{a v}$ subsides in favor of the $D$ eigenstate, and then, the latter will dictate the state vector output. This dynamical behavior is confirmed from (c) and (d), where under CW/CCW conditions (solid/dotted curves), switching to modes $\tilde{(+)}$ / $\tilde{(-)}$ $(D_{a v})$ and $\tilde{(-)}$ / $\tilde{(-)}$ $(D)$ takes place for T₁ and T₂, leading to ASC and SSC, respectively. (e)–(h) Trajectory parameters: ${\overset{̌}{g}}_{C} = 1.3, {\overset{̌}{σ}}_{C} = 0.1, {\overset{̌}{ε}}_{g} = 0.7, {\overset{̌}{ε}}_{σ} = 0.4, \bar{τ} = 1 / 64, κ T_{1} = 300, κ T_{2} = 3800$ ⁠. It can be deduced from (e) and (f) that $D_{a v} = \tilde{(+)}$ / $\tilde{(-)}$ and $D = \tilde{(-)}$ / $\tilde{(-)}$ ⁠, for CW/CCW evolution. Based then on (g) and (h), the open quantum system converts to the $D_{a v}$ and $D$ modes for periods T₁ and T₂, thereby resulting in ASC and SSC responses, correspondingly. (i)–(l) Trajectory parameters: ${\overset{̌}{g}}_{C} = 0.6, {\overset{̌}{σ}}_{C} = 0.2, {\overset{̌}{ε}}_{g} = {\overset{̌}{ε}}_{σ} = 0.4, \bar{τ} = - 1 / 4, κ T_{1} = 200, κ T_{2} = 2000$ ⁠. Panels (i) and (j) indicate that $D_{a v} = \tilde{(-)}$ and $D = \tilde{(+)}$ ⁠, for both encircling directions. Hence, switching to $\tilde{(-)}$ and $\tilde{(+)}$ mode branches should be observed for T₁ and T₂, respectively, which is verified from (k) and (l). Evidently, SSC occurs for the trajectory under study, with regard to both T₁ and T₂; however, additional numerical results for intermediate periods T suggest that ASC can still be attained [see the inset of (l) with κT₃ = 1000; the main graph and inset share the same vertical axis margins]. To interpret this, one should consider that a particular path, whether open or closed, generally yields distinct dynamics when traversed forward (CW) vs backward (CCW). As such, T_cr depends upon the path orientation, which implies the existence of intervals of T where evolution in the forward and backward directions does not entail conversion exclusively to the averagely $(D_{a v})$ or to the instantaneous $(D)$ dominant eigenstates. A specific example is demonstrated in the inset of (l): the dynamics at the end of the CW and CCW cycles are determined by the $D_{a v}$ $[\tilde{(-)}]$ and $D$ $[\tilde{(+)}]$ eigenmodes, respectively, therefore leading to asymmetric switching.

FIG. 5.

View large Download slide

Reconfigurable operation between the ASC and SSC switching regimes. The first two columns illustrate the closed trajectories $C^{R}$ within $R$ space (insets) along with the evolution of the eigenvalues ${\overset{̌}{λ}}^{\tilde{\pm}}$ in the complex domain (main panels), under clockwise (first column) and counterclockwise (second column) conditions. The last two columns depict the dynamics of the ratio factor for moderate (third column—T₁) and large (fourth column—T₂) values of the cycle period T (state conversion occurs in both occasions, but to modes $D_{a v}$ and $D$ ⁠, respectively, as we shall see). The colors and types of lines and symbols appearing in the various graphs have already been described in Fig. 4. (a)–(d) Trajectory parameters: ${\overset{̌}{g}}_{C} = 1, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = {\overset{̌}{ε}}_{σ} = 0.2, \bar{τ} = 1 / 16, κ T_{1} = 100, κ T_{2} = 2500$ ⁠. [Variable $\bar{τ} = τ / T$ denotes the CW $(\bar{τ} > 0)$ or CCW $(\bar{τ} < 0)$ normalized time shift from the left corner of $C^{R}$ ⁠.] For CW encircling [(a)], it can be inferred that $D_{a v} = \tilde{(+)}$ and $D = \tilde{(-)}$ ⁠, while in the CCW scenario [(b)], $D = D_{a v} = \tilde{(-)}$ ⁠. Consequently, and according to the analysis of Sec. III and Appendix B 3, conversion to the $D_{a v}$ mode is originally expected. Yet, as T is increased, there will be a critical transition point (T_cr) where $D_{a v}$ subsides in favor of the $D$ eigenstate, and then, the latter will dictate the state vector output. This dynamical behavior is confirmed from (c) and (d), where under CW/CCW conditions (solid/dotted curves), switching to modes $\tilde{(+)}$ / $\tilde{(-)}$ $(D_{a v})$ and $\tilde{(-)}$ / $\tilde{(-)}$ $(D)$ takes place for T₁ and T₂, leading to ASC and SSC, respectively. (e)–(h) Trajectory parameters: ${\overset{̌}{g}}_{C} = 1.3, {\overset{̌}{σ}}_{C} = 0.1, {\overset{̌}{ε}}_{g} = 0.7, {\overset{̌}{ε}}_{σ} = 0.4, \bar{τ} = 1 / 64, κ T_{1} = 300, κ T_{2} = 3800$ ⁠. It can be deduced from (e) and (f) that $D_{a v} = \tilde{(+)}$ / $\tilde{(-)}$ and $D = \tilde{(-)}$ / $\tilde{(-)}$ ⁠, for CW/CCW evolution. Based then on (g) and (h), the open quantum system converts to the $D_{a v}$ and $D$ modes for periods T₁ and T₂, thereby resulting in ASC and SSC responses, correspondingly. (i)–(l) Trajectory parameters: ${\overset{̌}{g}}_{C} = 0.6, {\overset{̌}{σ}}_{C} = 0.2, {\overset{̌}{ε}}_{g} = {\overset{̌}{ε}}_{σ} = 0.4, \bar{τ} = - 1 / 4, κ T_{1} = 200, κ T_{2} = 2000$ ⁠. Panels (i) and (j) indicate that $D_{a v} = \tilde{(-)}$ and $D = \tilde{(+)}$ ⁠, for both encircling directions. Hence, switching to $\tilde{(-)}$ and $\tilde{(+)}$ mode branches should be observed for T₁ and T₂, respectively, which is verified from (k) and (l). Evidently, SSC occurs for the trajectory under study, with regard to both T₁ and T₂; however, additional numerical results for intermediate periods T suggest that ASC can still be attained [see the inset of (l) with κT₃ = 1000; the main graph and inset share the same vertical axis margins]. To interpret this, one should consider that a particular path, whether open or closed, generally yields distinct dynamics when traversed forward (CW) vs backward (CCW). As such, T_cr depends upon the path orientation, which implies the existence of intervals of T where evolution in the forward and backward directions does not entail conversion exclusively to the averagely $(D_{a v})$ or to the instantaneous $(D)$ dominant eigenstates. A specific example is demonstrated in the inset of (l): the dynamics at the end of the CW and CCW cycles are determined by the $D_{a v}$ $[\tilde{(-)}]$ and $D$ $[\tilde{(+)}]$ eigenmodes, respectively, therefore leading to asymmetric switching.

FIG. 6.

View large Download slide

Slowly varying non-Hermitian dynamics when crossing special points in $R$ space, in terms of the limiting behavior of families of parametric trajectories. Our focus here lies on paths that include EP degeneracies, or whose endpoints satisfy the relation $R e ({\overset{̌}{λ}}^{\tilde{-}}) = R e ({\overset{̌}{λ}}^{\tilde{+}})$ (and thus mode $D$ is not well-defined). Through a number of illuminating examples, it will become apparent how to cleverly construct families of trajectories $(C^{R})$ ⁠, which gradually approach the path $(C^{R})$ whose adiabatic properties we wish to investigate (the $C^{R}$ curves illustrated in the various graphs are chosen not to be too close to $C^{R}$ ⁠, for visual purposes). The representation logic in the above diagrams follows the pattern of Figs. 4 and 5, with the difference being that solid/dotted lines will refer to $C^{R}$ / $C^{R}$ ⁠. Furthermore, the circle/triangle and square/diamond symbols in the ratio factor plots signify the eigenvector component ratio $r^{\tilde{+}}$ / $r^{\tilde{-}}$ at $\bar{t} = 1$ (t = T) for $C^{R}$ and $C^{R}$ ⁠, respectively (in all cases, either $r^{\tilde{+}}$ or $r^{\tilde{-}}$ is depicted depending on which mode eventually dominates). (a) and (b) Here, we are interested on the dynamics along SE $(C^{R})$ ⁠, where $R e ({\overset{̌}{λ}}^{\tilde{\pm}}) = 0$ at location E [see the main panel and the top inset of (a)]. Two different families of $C^{R}$ paths are examined: SK with K → E and MG with MG → SE. For SK and MG, it can be inferred from (a) that $R e ({\overset{̌}{λ}}^{\tilde{+}}) > R e ({\overset{̌}{λ}}^{\tilde{-}})$ at their respective end points (and at all points for that matter), and thus, $D = \tilde{(+)}$ ⁠. Moreover, $R e (⟨ {\overset{̌}{λ}}^{\tilde{+}} ⟩) > R e (⟨ {\overset{̌}{λ}}^{\tilde{-}} ⟩)$ resulting in $D_{a v} = \tilde{(+)}$ (gray outlined circles signifying $⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩$ are not shown, to avoid visual congestion). Consequently, the quantum state will most definitely convert to the $\tilde{(+)}$ mode for SK and MG. Taking the limits SK → SE (K → E) and MG → SE, the previous conclusion will also hold for SE, and this is verified from (b) [the faded r curve and hexagram symbols, indicating the value of $r^{\tilde{+}}$ at $\bar{t} = 1$ ⁠, both correspond to SK in panel (b)]. The coordinates of the various points in $R$ space and evolution period are given by S(1, 0.5), K(2/3, 1/6), E(0.5, 0), M(1, 0.7), G(0.5, 0.2), and κT = 60 (MG is shifted 0.2 units above SE). Regarding the bottom inset of panel (a), which depicts the eigenspectra plots for path SE in terms of the (±) notation, it has been included solely to facilitate the discussion in the caption of Fig. 3. (c) and (d) In the current simulation scenario, EP crossing is noticed for the main linear path $C^{R}$ ⁠. As such, an additional elliptic trajectory $C^{R}$ is illustrated in the inset of (c), which circumvents the EP, and therefore, the findings of Sec. III and Appendix B 3 should apply. Indeed, both $D$ and $D_{a v}$ are described by the $\tilde{(+)}$ branch, and thus, switching to the $\tilde{(+)}$ mode is expected for $C^{R}$ ⁠, as shown in (d). Considering the limit $C^{R} \to C^{R}$ ⁠, conversion to $\tilde{(+)}$ mode is also anticipated for $C^{R}$ ⁠, and this is confirmed from (d). The coordinates for S and E, evolution period, and elliptic curve parameters are provided by S(1.5, 0.5), E(0.5, −0.5), $κ T = 100, μ = π / 4, R_{g} = 0.5 \sqrt{2}, R_{σ} = 0.1 \sqrt{2}$ (μ, R_g, and R_σ are the rotation angle, major axis, and minor axis, respectively, of the half ellipse $C^{R}$ ⁠). (e)–(h) A more complex situation is examined here, where an EP is crossed not only during [(e) and (f)], but at the end of evolution as well [(g) and (h)]. Instead of studying the dynamical behavior along $C^{R}$ ⁠, we turn our attention to $C^{R}$ ⁠, which does not include any EPs [see the insets of (e) and (g)]. Emphasizing initially on (e) and (f), it can be deduced that $D_{a v} = \tilde{(+)}$ and $D = \tilde{(-)}$ at the end of the cycle. Hence, switching to mode $\tilde{(+)}$ is expected for moderate values of T [κT = 100, see the main panel of (f)], while in the superadiabatic case [κT = 2200, see the inset of (f)], the dynamics will be dictated by the eigenstate $\tilde{(-)}$ ⁠. Moving on to (g) and (h), we present solely the superadiabatic limit scenario, where mode $D$ prevails [only the gray solid outlined circles are shown in (g), as quantities $⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩$ are nearly identical for $C^{R}$ and $C^{R}$ ]. The instantaneous dominant eigenstate at point M is identified with mode $\tilde{(-)}$ ⁠, which will describe the quantum state response [see (h), with κT = 500]. All the aforementioned observations can now be generalized to $C^{R}$ after taking the limit $C^{R} \to C^{R}$ ⁠, and this is verified from the solid curves in (f) and (h). The parameters employed in (e) and (f) $(\bar{τ} = - 1 / 8)$ and (g) and (h) $(\bar{τ} = - 1 / 2)$ are as follows for $C^{R}$ $(C^{R})$ ⁠: ${\overset{̌}{g}}_{C} = 0.9, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = 0.1$ (0.3), ${\overset{̌}{ε}}_{σ} = 0.1$ ⁠.

FIG. 6.

View large Download slide

Slowly varying non-Hermitian dynamics when crossing special points in $R$ space, in terms of the limiting behavior of families of parametric trajectories. Our focus here lies on paths that include EP degeneracies, or whose endpoints satisfy the relation $R e ({\overset{̌}{λ}}^{\tilde{-}}) = R e ({\overset{̌}{λ}}^{\tilde{+}})$ (and thus mode $D$ is not well-defined). Through a number of illuminating examples, it will become apparent how to cleverly construct families of trajectories $(C^{R})$ ⁠, which gradually approach the path $(C^{R})$ whose adiabatic properties we wish to investigate (the $C^{R}$ curves illustrated in the various graphs are chosen not to be too close to $C^{R}$ ⁠, for visual purposes). The representation logic in the above diagrams follows the pattern of Figs. 4 and 5, with the difference being that solid/dotted lines will refer to $C^{R}$ / $C^{R}$ ⁠. Furthermore, the circle/triangle and square/diamond symbols in the ratio factor plots signify the eigenvector component ratio $r^{\tilde{+}}$ / $r^{\tilde{-}}$ at $\bar{t} = 1$ (t = T) for $C^{R}$ and $C^{R}$ ⁠, respectively (in all cases, either $r^{\tilde{+}}$ or $r^{\tilde{-}}$ is depicted depending on which mode eventually dominates). (a) and (b) Here, we are interested on the dynamics along SE $(C^{R})$ ⁠, where $R e ({\overset{̌}{λ}}^{\tilde{\pm}}) = 0$ at location E [see the main panel and the top inset of (a)]. Two different families of $C^{R}$ paths are examined: SK with K → E and MG with MG → SE. For SK and MG, it can be inferred from (a) that $R e ({\overset{̌}{λ}}^{\tilde{+}}) > R e ({\overset{̌}{λ}}^{\tilde{-}})$ at their respective end points (and at all points for that matter), and thus, $D = \tilde{(+)}$ ⁠. Moreover, $R e (⟨ {\overset{̌}{λ}}^{\tilde{+}} ⟩) > R e (⟨ {\overset{̌}{λ}}^{\tilde{-}} ⟩)$ resulting in $D_{a v} = \tilde{(+)}$ (gray outlined circles signifying $⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩$ are not shown, to avoid visual congestion). Consequently, the quantum state will most definitely convert to the $\tilde{(+)}$ mode for SK and MG. Taking the limits SK → SE (K → E) and MG → SE, the previous conclusion will also hold for SE, and this is verified from (b) [the faded r curve and hexagram symbols, indicating the value of $r^{\tilde{+}}$ at $\bar{t} = 1$ ⁠, both correspond to SK in panel (b)]. The coordinates of the various points in $R$ space and evolution period are given by S(1, 0.5), K(2/3, 1/6), E(0.5, 0), M(1, 0.7), G(0.5, 0.2), and κT = 60 (MG is shifted 0.2 units above SE). Regarding the bottom inset of panel (a), which depicts the eigenspectra plots for path SE in terms of the (±) notation, it has been included solely to facilitate the discussion in the caption of Fig. 3. (c) and (d) In the current simulation scenario, EP crossing is noticed for the main linear path $C^{R}$ ⁠. As such, an additional elliptic trajectory $C^{R}$ is illustrated in the inset of (c), which circumvents the EP, and therefore, the findings of Sec. III and Appendix B 3 should apply. Indeed, both $D$ and $D_{a v}$ are described by the $\tilde{(+)}$ branch, and thus, switching to the $\tilde{(+)}$ mode is expected for $C^{R}$ ⁠, as shown in (d). Considering the limit $C^{R} \to C^{R}$ ⁠, conversion to $\tilde{(+)}$ mode is also anticipated for $C^{R}$ ⁠, and this is confirmed from (d). The coordinates for S and E, evolution period, and elliptic curve parameters are provided by S(1.5, 0.5), E(0.5, −0.5), $κ T = 100, μ = π / 4, R_{g} = 0.5 \sqrt{2}, R_{σ} = 0.1 \sqrt{2}$ (μ, R_g, and R_σ are the rotation angle, major axis, and minor axis, respectively, of the half ellipse $C^{R}$ ⁠). (e)–(h) A more complex situation is examined here, where an EP is crossed not only during [(e) and (f)], but at the end of evolution as well [(g) and (h)]. Instead of studying the dynamical behavior along $C^{R}$ ⁠, we turn our attention to $C^{R}$ ⁠, which does not include any EPs [see the insets of (e) and (g)]. Emphasizing initially on (e) and (f), it can be deduced that $D_{a v} = \tilde{(+)}$ and $D = \tilde{(-)}$ at the end of the cycle. Hence, switching to mode $\tilde{(+)}$ is expected for moderate values of T [κT = 100, see the main panel of (f)], while in the superadiabatic case [κT = 2200, see the inset of (f)], the dynamics will be dictated by the eigenstate $\tilde{(-)}$ ⁠. Moving on to (g) and (h), we present solely the superadiabatic limit scenario, where mode $D$ prevails [only the gray solid outlined circles are shown in (g), as quantities $⟨ {\overset{̌}{λ}}^{\tilde{\pm}} ⟩$ are nearly identical for $C^{R}$ and $C^{R}$ ]. The instantaneous dominant eigenstate at point M is identified with mode $\tilde{(-)}$ ⁠, which will describe the quantum state response [see (h), with κT = 500]. All the aforementioned observations can now be generalized to $C^{R}$ after taking the limit $C^{R} \to C^{R}$ ⁠, and this is verified from the solid curves in (f) and (h). The parameters employed in (e) and (f) $(\bar{τ} = - 1 / 8)$ and (g) and (h) $(\bar{τ} = - 1 / 2)$ are as follows for $C^{R}$ $(C^{R})$ ⁠: ${\overset{̌}{g}}_{C} = 0.9, {\overset{̌}{σ}}_{C} = 0, {\overset{̌}{ε}}_{g} = 0.1$ (0.3), ${\overset{̌}{ε}}_{σ} = 0.1$ ⁠.

We shall turn now our attention to the simulation results depicted in the third column of Fig. 5, where conversion to the averagely dominant mode $D_{a v}$ occurs. [We need not further analyze the graphs of Fig. 4 with respect to $D_{a v}$ ⁠, as either superadiabaticity is reached and thus it cannot emerge (second and third rows of Fig. 4), or simply $D = D_{a v}$ (first row of Fig. 4) and therefore $D_{a v}$ does not alter in any way the resulting dynamics.] Relative to the fourth column, period T has been reduced, so as the topology of the entire path $C^{R}$ has an effect on the identity of the dominating eigenstate. Consequently, instead of getting a SSC-type response for all parameter configurations, asymmetric switching is noticed in the first and second rows of Fig. 5. Notably, conversion to the $D_{a v}$ mode can still manifest in the presence of EPs, as evidenced in Figs. 6(c)–6(f). {Figures 6(e) and 6(f) [see the main graph of Fig. 6(f)] provide a more illustrative example of the latter, since then $D \neq D_{a v}$ ⁠.}

Bringing it altogether, an important and previously overlooked dynamic trait of slow non-Hermitian cycling is unveiled: the capacity to promote reconfigurable ASC/SSC operation by tuning the overall period of evolution. For instance, if T is increased for the scenarios demonstrated in the first two rows of Fig. 5, the switching response shifts from ASC (third column) to SSC (fourth column). An even more intricate switching pattern is attained in the last row of the same figure, as T grows larger: SSC [Fig. 5(k)] → ASC [the inset of Fig. 5(l)] → SSC [the main panel of Fig. 5(l)]. To interpret these peculiar (yet extremely useful) dynamical features, we need to consider that traversing a path $C^{R}$ in the forward (CW) and backward (CCW) manner will produce different state vector evolutions [this can be directly seen via Eq. (54), where in general Ξ_NΞ_N−1⋯Ξ₀ ≠ Ξ₀⋯Ξ_N−1Ξ_N given the non-commutativity of matrix multiplication], except from very special cases [e.g., if H(t) = H(T − t), M(t) = M(T − t) in Eqs. (1) and (2), or equivalently if Ξ_n = Ξ_N−n, H_n = H_N−n, M_n = M_N−n in Eq. (54), or when ASC occurs in cyclic non-Hermitian settings]. In this regard, and assuming that conditions for the emergence of $D_{a v}$ are met, $C^{R}$ will be characterized by two distinct critical values, $T_{c r_{1}}$ and $T_{c r_{2}}$ ⁠, signifying the thresholds below/above which the quantum state converts to $D_{a v} / D$ ⁠, with respect to the two opposite traversal directions. The preceding suggests that a maximum of three different (and neighboring) ranges of T can arise, each resulting in ASC or SSC responses. The exact number of such intervals varies according to whether $D \neq D_{a v}$ or $D = D_{a v}$ under CW/CCW conditions, and if mode $D_{a v}$ is supported or not. Along these lines, a variety of switching patterns can be acquired, depending on the chosen topology of $C^{R}$ ⁠.

The foregoing arguments ultimately lead to the proposition for a reconfigurable optical switch/omnipolarizer based on fiber loop networks. Such synthetic photonic lattices, which were outlined in Sec. II B, allow for the realization of the continuous two-level Schrödinger equation depicted in Eq. (1) or (2), but in a discrete time fashion [see Eq. (54)]. (In what follows, we shall briefly describe the operation of fiber loop setups, but for more technical details, the reader should resort to the supplementary material of Ref. 22 or Ref. 28.) The state vector will correspond to the polarization components E_x, E_y of the propagating field, which can be coupled via a polarization controller (implemented through a variable retardation wave-plate or a birefringent fiber) and individually manipulated by intensity (amplifying/attenuating E_x, E_y) and phase modulators (imparting phase on E_x, E_y). Beam splitters separate the polarization components for independent modulation, while beam combiners merge them back together. Each round trip along the fiber loop marks an incremental time step (t_n → t_n+1, t_n+1 − t_n = Δ), and thus, the adiabatic rate will be dictated by the overall number of round trips N (β = 1/T ∝ 1/N, where T = NΔ and Δ is a fixed small time step). The great degree of tunability inherent in these discrete optical configurations enables the implementation of diverse $C^{R}$ topologies, thereby giving rise to a wide range of switching patterns [the exact type of response (ASC or SSC) can be specified, after appropriately adjusting N]. While this paragraph focuses on a particular class of photonic networks, our findings have broader applicability to quantum or other physical settings as long as experimental platforms with comparable levels of modulation flexibility are available.

To summarize, we would like to make the following remarks:

The adiabatic transport features underpinning the non-Hermitian evolution can be determined by examining the topology of parametric contour $C^{R}$ ⁠, either at (or near) the end, or in an average fashion. In the former (latter) occasion, which arises for sufficiently large (moderate) values of T, we need to compare quantities $R e [λ^{\tilde{\pm}} (T)]$ $[R e ({⟨ λ^{\tilde{\pm}} ⟩}_{T})]$ to identify the prevalent eigenstate. The above criteria pertain to both open and closed trajectories in $R$ space and are consistent with our recent study in Ref. 22, where the state vector dynamics of discrete non-Hermitian cyclic arrangements were analyzed but only in the superadiabatic limit.
Exceptional point encirclement is neither a necessary [see Figs. 4(e) and 4(f)] nor a sufficient [see Figs. 4(c) and 4(d)] condition for asymmetric state conversion to appear. This holds true for symmetric switching as well [compare Figs. 4(a) and 4(b) with Figs. 4(c) and 4(d)]. In Appendix C, a separate investigation has been conducted as to the effect of EP encirclement on the resulting switching performance.
The conclusions of the current analysis will continue to apply in the presence of points where eigenstates $D$ or $D_{a v}$ are not well-defined, as inferred from Fig. 6. The identities of $D$ and $D_{a v}$ can then be specified, by inspecting quantities $R e [λ^{\tilde{\pm}} (t)]$ and $R e ({⟨ λ^{\tilde{\pm}} ⟩}_{t})$ in the limit t → T. Furthermore, the crossing of EPs is not prohibitive for any state conversion effects to take place (of course, operation at such degeneracies is highly unlikely in actual experimental situations, owing to their inherent sensitivity^34–36). All the preceding circumstances are reviewed in a more detailed manner within Appendix B 3, in terms of the mode population equations.
The technique, presented in Fig. 6, can be utilized to investigate/verify the dynamics along arbitrary parametric paths, at least when non-Hermitian conditions are considered and mode switching is noticed (this in turn signifies the robustness of non-Hermitian switching processes). Yet, if mode switching does not occur, the output state will typically vary significantly even with small changes in $C^{R}$ [the previous can be deduced, for instance, from Figs. 3(d) and 3(e), where the state vector fluctuates considerably with time], and then, the methodology of Fig. 6 may not offer distinct advantages.
In the current section, we explored the potential of attaining reconfigurable ASC/SSC operation within cyclic non-Hermitian environments, by exploiting the adiabatic state transition to modes $D$ or $D_{a v}$ contingent upon the magnitude of T (this subsequently led to our proposal for reconfigurable omnipolarizer devices). An analogous functionality is supported in acyclic non-conservative settings, but then, it is not pertinent to discuss about (a)symmetric state conversion, as the starting and ending points in $R$ space do not overlap. Essentially, in the case of open parametric trajectories, we examine the state vector response as they are traversed in the forward and backward directions, which can also yield a rich variety of switching motifs.
The $P T$ -symmetric-like Hamiltonian, introduced in Sec. II B and used in the numerical (and some of the theoretical) computations, serves as a prototypical model upon which much of the existing literature on non-Hermitian structures is based. However, the implications of the present study extend beyond such a framework, to non-Hermitian Hamiltonians of arbitrary functional forms, as evidenced from our systematic analytical derivations. Furthermore, the reported switching to the $D$ or $D_{a v}$ modes is independent of the initial excitation (here, we assumed random initial conditions, yet the results remain unaffected if the instantaneous eigenstates get excited) and the precise shape of the traversed path in $R$ space [in the various numerical scenarios, we emphasized on linear parametric trajectories and linear time dependencies, but similar observations hold for alternative path profiles or temporal evolutions (see Fig. 1)].

V. CONCLUSIONS AND OUTLOOK

In this paper, the dynamics of time-dependent two-level open quantum systems, both cyclic and acyclic, have been rigorously investigated in the slowly varying limit for parametric trajectories of arbitrary geometries. Such an analysis has been conducted both in a discrete (state vector component ratio) and in a continuous time fashion (mode population equations), and it has been shown that conversion to the least dissipative instantaneous eigenstate $(D)$ is expected for superadiabatic processes (special emphasis is then given to the topology of the respective eigenmode branches toward the end of the parametric evolution). This has been confirmed in the particular case of traceless and symmetric non-Hermitian (or $P T$ -symmetric-like) Hamiltonian configurations, where exact solutions have been retrieved for two distinct types of time dependencies (linear/hyperbolic—Airy/modified Bessel solutions). Our asymptotic calculations for very large evolution periods have clearly demonstrated the role of the Stokes phenomenon in the observed mode switching behavior. While a two-dimensional $R$ space was assumed in the aforementioned scenario (defined by the gain and detuning coefficients), the main analysis and conclusions of our work (Secs. II A, III, and IV) are not constrained by the dimensionality of the parameter space.

An important outcome of the current study is that there exists a transitional range of values for the adiabatic rate β, below and above which conversion still takes place but not necessarily to the same eigenstate. This led to the introduction of the averagely dominant eigenstate $(D_{a v})$ ⁠, whose identity is determined after examining the topology of the eigenmode branches as a whole and not only toward their ending points (the competition for dominance between $D$ and $D_{a v}$ has been illustrated within Appendix B 3 via a detailed inspection of the mode population equations, while also highlighting the significance of the non-adiabatic mode couplings). All our findings have been eventually applied in the case of closed trajectories in $R$ space, to provide a reasoning behind the recently reported SSC and ASC effects in cyclic non-conservative Hamiltonian arrangements. Moreover, we have proposed a reconfigurable optical switch/omnipolarizer, which under cyclic (acyclic) conditions can exhibit either SSC or ASC behaviors (can switch either to the $D$ or to the $D_{a v}$ modes) by tuning the overall period of evolution.

To conclude this discussion, even though the great majority of the aspects concerning the adiabatic properties of two-mode non-Hermitian Hamiltonian settings have been elucidated, there are still certain topics that need to be separately addressed. One of them is to derive more accurate estimates for the critical threshold of β that signifies (i) the transition of the quantum state response from $D_{a v}$ to $D$ (β_cr) and (ii) the onset of any state conversion phenomena (β_ad). Along these lines, a question arises as to whether it is possible under certain circumstances for β_ad < β_cr, which consequently implies that $D_{a v}$ cannot emerge (the quantum state can convert then only to the $D$ eigenstate). Furthermore, it should be noted that while the physical intuition and results presented here can be readily extended to higher (and even infinite)-dimensional Hilbert spaces, the underlying mathematical techniques should be appropriately adapted in such situations. Actually, this will be the subject of forthcoming investigations; already though, some characteristic examples involving three-level non-Hermitian Hamiltonians have been provided in the supplementary material. Finally, of great practical interest is to explore the impact of any nonlinear (or other external perturbation) effects on the resulting adiabatic transport dynamics, as this will enable a more systematic evaluation of the performance and robustness of mode switching platforms inspired from non-Hermitian quantum theory.

SUPPLEMENTARY MATERIAL

The supplementary material provides details as to how the developed methodologies can be extended to higher-dimensional Hamiltonian settings. In this respect, the specific case of 3 × 3 non-Hermitian Hamiltonians has been analytically examined in the superadiabatic limit, in terms of both the mode population and the state vector component ratio dynamics (in the latter instance, a two-dimensional generalization of the G-transform has been utilized). Switching to the instantaneous dominant eigenmode D has been shown then, under both cyclic and acyclic conditions, given an imbalance among the real parts of the system eigenvalues λ(t).

ACKNOWLEDGMENTS

N.S.N. was supported by the Bodossaki Foundation.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Nicholas S. Nye: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Nikolaos V. Kantartzis: Validation (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available within the article.

APPENDIX A: SLOWLY VARYING BEHAVIOR OF GENERALIZED QUANTUM HAMILTONIAN SETTINGS IN TERMS OF THE WKB METHOD

1. Calculating the eigenvalues and the left and right eigenvectors of a two-level Hamiltonian system

The left and right eigenvalue problems for the system matrix M(t) can be written, respectively, as

⟨v_{L} (t)| [λ (t) I - M (t)] = 0,

(A1a)

[λ (t) I - M (t)] |v_{R} (t)⟩ = 0 .

(A1b)

[Within the main text, when we refer to eigenvectors, we typically imply the right eigenvectors

|v_{R} (t)⟩

⁠.] The associated characteristic equation for M(t) becomes

\det [λ (t) I - M (t)] = 0 \Leftrightarrow λ^{2} (t) - t r [M (t)] λ (t) + \det [M (t)] = 0,

(A2)

which accepts as solutions the eigenvalues

λ^{\pm} (t) = \frac{t r [M (t)] \pm i \sqrt{4 \det [M (t)] - {t r}^{2} [M (t)]}}{2}

(A3)

and is also described by the subsequent Vieta’s formulas,

λ^{+} (t) + λ^{-} (t) = t r [M (t)],

(A4a)

λ^{+} (t) λ^{-} (t) = \det [M (t)] .

(A4b)

From Eq. , it can be seen that the ratio of the components of the left

⟨v_{L}^{\pm} (t)| = [v_{L, x}^{\pm} (t) v_{L, y}^{\pm} (t)]

and right

|v_{R}^{\pm} (t)⟩ = {[v_{R, x}^{\pm} (t) v_{R, y}^{\pm} (t)]}^{T}

eigenvectors will be given as

r_{L}^{\pm} (t) = \frac{v_{L, y}^{\pm} (t)}{v_{L, x}^{\pm} (t)} = \frac{λ^{\pm} (t) - m^{11} (t)}{m^{21} (t)} = \frac{m^{12} (t)}{λ^{\pm} (t) - m^{22} (t)},

(A5a)

r_{R}^{\pm} (t) = \frac{v_{R, y}^{\pm} (t)}{v_{R, x}^{\pm} (t)} = \frac{λ^{\pm} (t) - m^{11} (t)}{m^{12} (t)} = \frac{m^{21} (t)}{λ^{\pm} (t) - m^{22} (t)} .

(A5b)

Note that in the forthcoming derivations, we shall assign

r_{L}^{\pm} (t) = \frac{m^{12} (t)}{m^{21} (t)} r^{\pm} (t), r_{R}^{\pm} (t) = r^{\pm} (t),

(A6)

but in the present subappendix, we will keep using variables

r_{L}^{\pm} (t), r_{R}^{\pm} (t)

⁠. Following now Eq. , we have that

⟨v_{L}^{\pm} (t)| = [1 r_{L}^{\pm} (t)], |v_{R}^{\pm} (t)⟩ = {[1 r_{R}^{\pm} (t)]}^{T},

(A7)

which lead to the subsequent inner products,

⟨v_{L}^{\pm} (t)| v_{R}^{\pm} (t)⟩ = 1 + r_{L}^{\pm} (t) r_{R}^{\pm} (t),

(A8a)

\begin{align} ⟨v_{L}^{\mp} (t)| v_{R}^{\pm} (t)⟩ & = 1 + r_{L}^{\mp} (t) r_{R}^{\pm} (t) = 1 + \frac{λ^{\mp} (t) - m^{11} (t)}{λ^{\pm} (t) - m^{22} (t)} \\ = \frac{λ^{\pm} (t) + λ^{\mp} (t) - m^{11} (t) - m^{22} (t)}{λ^{\pm} (t) - m^{22} (t)} = 0 . \end{align}

(A8b)

As the above relations reveal, the left and right eigenvectors of M(t) form a biorthogonal basis. A biorthonormal basis can also be acquired as follows:

⟨V_{L}^{\pm} (t)| = \frac{⟨v_{L}^{\pm} (t)|}{1 + r_{L}^{\pm} (t) r_{R}^{\pm} (t)}, |V_{R}^{\pm} (t)⟩ = |v_{R}^{\pm} (t)⟩,

(A9a)

⟨V_{L}^{\pm} (t)| V_{R}^{\pm} (t)⟩ = 1, ⟨V_{L}^{\mp} (t)| V_{R}^{\pm} (t)⟩ = 0,

(A9b)

and will be frequently employed in our analysis. It should be highlighted that the choice of the normalization scheme just introduces a gauge degree of freedom, which does not have an effect in any measurable (gauge-invariant) quantities [this will be further discussed in Appendix B 2 via the particular example of (a)cyclic geometric phases and their effect on the resulting state vector].

Of interest is to also provide the conditions under which spectrum degeneracies arise in two-level Hamiltonian systems. These are coined as either diabolic (DPs) or exceptional (EPs) points in the Hermitian and non-Hermitian cases, respectively, but in the current investigation, we are particularly interested in the latter. According then to Eq. ,

λ_{EP}^{+} = λ_{EP}^{-}

can apply as long as

{tr}^{2} (M_{EP}) = 4 \det (M_{EP}) \Leftrightarrow m_{EP}^{11} - m_{EP}^{22} = \pm 2 i \sqrt{m_{EP}^{12} m_{EP}^{21}} .

(A10)

(The ± sign in the previous formula should not be confused with earlier equations, where it signified the two distinct eigenstates of a two-mode Hamiltonian setting. Moreover, we have omitted the dependence on time t and instead introduced the subscript EP, to denote operation at EPs). The degenerate eigenvalue λ_EP

(= λ_{EP}^{\pm})

will be given by

λ_{EP} = \frac{tr (M_{EP})}{2} = \frac{m_{EP}^{11} + m_{EP}^{22}}{2},

(A11)

and the respective left and right eigenvector component ratios from

r_{L}^{EP} = \mp i \frac{\sqrt{m_{EP}^{12} m_{EP}^{21}}}{m_{EP}^{21}}, r_{R}^{EP} = \mp i \frac{\sqrt{m_{EP}^{12} m_{EP}^{21}}}{m_{EP}^{12}} .

(A12)

[We have been careful in the above expression not to simply write

r_{L}^{EP} = \mp i \sqrt{m_{EP}^{12} / m_{EP}^{21}}

and

r_{R}^{EP} = \mp i \sqrt{m_{EP}^{21} / m_{EP}^{12}}

⁠, as this requires relation

\sqrt{z_{1} z_{2}} = \sqrt{z_{1}} \sqrt{z_{2}}

with z_1,2 arbitrary complex variables to uphold, but this is not always true when adhering to the principal (or any other branch) of the square root function (see Appendix A 5 and Sec. II B for more extensive calculations with principal square roots).] It can be attained now that

⟨v_{L}^{EP}| v_{R}^{EP}⟩ = 1 + r_{L}^{EP} r_{R}^{EP} = 0,

(A13)

which can alternatively be inferred from Eq. after substituting

⟨v_{L}^{\mp} (t)| \to ⟨v_{L}^{EP}|

and

|v_{R}^{\pm} (t)⟩ \to |v_{R}^{EP}⟩

(left and right eigenvectors coalesce). From Eqs. and , it becomes evident that eigenvector normalization cannot take place at the location of the EP degeneracies.

Finally, we shall repeat the most essential of the above formulas in the special case of traceless and symmetric non-Hermitian structures [m¹¹(t) = −m²²(t) = f(t), m¹²(t) = m²¹(t) = iκ, κ: constant and positive coupling strength], which were introduced in Sec. II B. The eigenvalues will then be given by

λ^{\pm} (t) = \pm i \sqrt{κ^{2} - f^{2} (t)},

(A14)

while the eigenvector component ratios for the left and right eigenvalue problems become both equal to

r_{L}^{\pm} (t) = r_{R}^{\pm} (t) = \pm \sqrt{1 - \frac{f^{2} (t)}{κ^{2}}} + i \frac{f (t)}{κ},

(A15)

and we have also used that

κ \in R^{+}

[if κ > 0, then

\sqrt{κ^{2} - f^{2} (t)} / κ = \sqrt{1 - f^{2} (t) / κ^{2}}

]. At the location of the EP degeneracies, it will be true that

λ_{EP} = 0 \Leftrightarrow f_{EP} = \mp κ .

(A16)

Furthermore, the following relation will hold:

r_{L}^{EP} = r_{R}^{EP} = \mp i .

(A17)

We shall conclude the present section by making the following comments:

Here and overall in our analysis, we emphasize on the system matrix M(t) instead of the Hamiltonian matrix H(t), owing to the non-Hermiticity of the quantum arrangements under study [H^†(t) ≠ H(t) or equivalently M^†(t) ≠ −M(t), given that M(t) = −iH(t)]. This allows us to turn our attention to the real part of the eigenvalues λ^±(t) of M(t), instead of the imaginary part of the eigenvalues $λ_{H}^{\pm} (t)$ of H(t) (known in the literature as eigenenergies), and thus avoid the multiple appearances of the imaginary unit i, which can be rather confusing. Of course, since M(t) and H(t) are linearly related [M(t) = −iH(t)], the same will be true for the associated eigenvalues $[λ^{\pm} (t) = - i λ_{H}^{\pm} (t)]$ with the left and right eigenvectors remaining unchanged. (Similarly, for the discrete evolution studied in Sec. III and further analyzed in Appendixes D and E, it applies that Ξ_n = I + M_nΔ = I − iH_nΔ, which in turn implies that matrices Ξ_n, M_n, H_n share the same eigenvectors and $λ_{Ξ_{n}}^{\pm} = 1 + λ_{n}^{\pm} Δ = 1 - i λ_{H_{n}}^{\pm} Δ$ for the respective eigenvalues.)
In Hermitian settings [H^†(t) = H(t) or M^†(t) = −M(t)], the eigenvalues λ^±(t) become imaginary [eigenenergies $λ_{H}^{\pm} (t)$ are real], while the left $[⟨v_{L}^{\pm} (t)|]$ and right $[|v_{R}^{\pm} (t)⟩]$ eigenvectors are conjugate transposes of each other,
$\begin{gathered} ⟨v_{L}^{\pm} (t)| M (t) = ⟨v_{L}^{\pm} (t)| λ^{\pm} (t) \Leftrightarrow \\ M^{†} (t) |v_{L}^{\pm} (t)⟩ = {[λ^{\pm} (t)]}^{*} |v_{L}^{\pm} (t)⟩ \Leftrightarrow \\ M (t) |v_{L}^{\pm} (t)⟩ = - {[λ^{\pm} (t)]}^{*} |v_{L}^{\pm} (t)⟩ \Leftrightarrow \\ {[λ^{\pm} (t)]}^{*} = - λ^{\pm} (t), |v_{L}^{\pm} (t)⟩ = |v_{R}^{\pm} (t)⟩ . \end{gathered}$
(A18)
{In the above formula, ${[λ^{\pm} (t)]}^{*}$ stands for the complex conjugate of λ^±(t).} Moreover, for spectrum degeneracies or diabolic points to arise $(λ_{DP}^{+} = λ_{DP}^{-})$ ⁠, it comes from Eq. (A10) that
$\begin{gathered} m_{DP}^{11} - m_{DP}^{22} = \pm 2 i \sqrt{m_{DP}^{12} m_{DP}^{21}} \Leftrightarrow \\ i (p - q) = \pm 2 i \sqrt{- | m_{DP}^{12} |^{2}} \Leftrightarrow \\ p - q = \pm 2 i | m_{DP}^{12} |, \end{gathered}$
(A19)
where we have applied relations $m_{DP}^{11} = - {(m_{DP}^{11})}^{*} = i p, m_{DP}^{22} = - {(m_{DP}^{11})}^{*} = i q$ ⁠, with $p, q \in R$ and $m_{DP}^{21} = - {(m_{DP}^{21})}^{*}$ (⁠ $M_{DP}^{†} = - M_{DP}$ ⁠, while index DP signifies operation at the diabolic point). Equation (A19) can hold iff p = q and $m_{DP}^{12} = m_{DP}^{21} = 0$ or equivalently M_DP = ipI (and H_DP = −pI for the Hamiltonian matrix). In this case, the system eigenvalues coincide at λ_DP = ip, but a fundamental difference can be noticed with respect to non-Hermitian degeneracies: whereas at EPs, a coalescence of eigenvectors takes place, this does not occur at DPs. Stated in terms of vector spaces, the dimensionality of the eigenspace (i.e., the vector space formed by the eigenvectors associated with a particular eigenvalue) at EPs is one, while at DPs is two, at least for the two-mode systems that we examine. [The only circumstance in which the dimensionality of the eigenspace at a non-Hermitian degeneracy is two (deeming it thus as a DP with non-coalescent eigenvectors) is when the system matrix M/Hamiltonian H becomes proportional to the identity matrix by some real/imaginary or in general complex quantity—such a marginal scenario should not concern us, as it has no effect on the ratio of the state vector components (which is of interest for us), and this can be directly deduced via the state vector discrete evolution equations depicted in Eq. (54).]
To more efficiently study the dynamical behavior of quantum arrangements in the quasistatic limit (T → ∞, T: period of evolution), the various system parameters (e.g., gain and detuning coefficients for the $P T$ -symmetric-like Hamiltonian of Sec. II B) are typically written directly as functions of the scaled time coordinate $\bar{t} = t / T$ (this entails that if t and T increase multiplicatively by the same factor, then the system parameters remain unaltered). This in turn implies that matrices H, M will depend directly on $\bar{t}$ ⁠, and the same will also be true for their respective eigenvalues and eigenvectors (the state vector will of course depend on both $\bar{t}, T$ ⁠, and thus, the aforementioned multiplicative rule will not hold then, as the subsequent WKB analysis will show).
Pertaining to $P T$ -symmetric-like settings, at many instances (particularly in numerical computations), a normalization scheme with respect to the coupling constant κ is employed. This arises after recasting the Hamiltonian evolution formula in Eq. (14) of Sec. II B as
$\frac{d |u (t)⟩}{d (κ t)} = [\begin{matrix} f (t) / κ & i \\ i & - f (t) / κ \end{matrix}] |u (t)⟩ .$
(A20)
The state vector dynamics can then be studied in terms of variables $\overset{̌}{t} = κ t, \overset{̌}{f} = f / κ$ $(\overset{̌}{T} = κ T, \overset{̌}{Δ} = κ Δ, \overset{̌}{g} = g / κ, \overset{̌}{σ} = σ / κ, {\overset{̌}{ε}}_{g, σ} = ε_{g, σ} / κ)$ ⁠. The state vector itself, its components, and their ratio remain unaffected (the only change is that now they are written as functions of $\overset{̌}{t}$ rather than t). The new system and Hamiltonian matrices read as $\overset{̌}{M} = M / κ, \overset{̌}{H} = H / κ$ ⁠, which implies for the eigenvalues that ${\overset{̌}{λ}}^{\pm} = λ^{\pm} / κ$ ⁠. Moreover, regarding the previous bullet point, it applies that $\bar{t} = \overset{̌}{t} / \overset{̌}{T} = t / T$ ⁠, and thus, the scaled time coordinate will not vary if we first normalize with respect to κ.

2. Retrieving the dynamical, geometric, and higher-order phase factors via WKB theory

We shall now put the WKB method into practice and study the behavior of the differential equation that the state vector component a(t) satisfies [see Eq. of Sec. II A] in the asymptotic limit T → ∞ [analogous conclusions can be drawn if instead we emphasized on component b(t)]. In this vein, the scaled temporal coordinate

\bar{t} = t / T

(t \in [0, T], \bar{t} \in [0, 1])

shall be employed and Eq. can then be reexpressed as

\begin{align} \frac{1}{T^{2}} \frac{d^{2} a^{μ} (\bar{t})}{d {\bar{t}}^{2}} - \{\frac{1}{T} t r [M (\bar{t})] + \frac{1}{T^{2}} / \frac{d m^{12} (\bar{t})}{d \bar{t}} m^{12} (\bar{t})\} \frac{d a^{μ} (\bar{t})}{d \bar{t}} \\ + \{\det [M (\bar{t})] + \frac{1}{T} m^{11} (\bar{t}) / \frac{d m^{12} (\bar{t})}{d \bar{t}} m^{12} (\bar{t}) \\ - \frac{1}{T} \frac{d m^{11} (\bar{t})}{d \bar{t}}\} a^{μ} (\bar{t}) = 0 . \end{align}

(A21)

The WKB ansatz assumes the following form:

a^{μ} (\bar{t}) \sim e^{μ \sum_{j = 0}^{\infty} S_{j} (\bar{t}) / μ^{j}} = e^{μ S_{0} (\bar{t}) + S_{1} (\bar{t}) + S_{2} (\bar{t}) / μ + O (μ^{- 2})},

(A22)

where

μ \in R^{+}

goes to infinity (the asymptotic scaling of μ with respect to T will be subsequently determined) and

a^{μ} (\bar{t})

signifies a dependence on both μ and

\bar{t}

[as opposed to the elements of matrix

M (\bar{t})

and terms

S_{j} (\bar{t})

⁠, which depend solely on

\bar{t}

]. The asymptotic notation ∼ in Eq. essentially signifies the divergent nature of the WKB approximation as μ → ∞, owing to the presence of the term

μ S_{0} (\bar{t})

[unless

S_{0} (\bar{t}) \equiv 0

]—this comes as no surprise, since for increasingly large T, the order of the differential equation in Eq. progressively reduces and as such the solution

a^{μ} (\bar{t})

exhibits an abrupt behavior in the asymptotic limit.

Given the aforementioned discussion, we require for

\sum S_{j} (\bar{t}) / μ^{j - 1}

to be an asymptotic (instead of strictly convergent) series in μ as μ → ∞ uniformly for all

\bar{t}

⁠. From d’Alembert’s ratio test, it should be true then that

| S_{j + 1} (\bar{t}) / S_{j} (\bar{t}) | ≪ μ

⁠, or equivalently that

S_{j + 1} (\bar{t}) / S_{j} (\bar{t})

is a bounded function of

\bar{t}

⁠. In typical scenarios, this might apply up to some value of j (let us say

J

⁠), above which

S_{j} (\bar{t}) / μ^{j - 1}

starts to increase in absolute value.³⁷ Hence, the optimal rule suggests truncating the WKB series at

S_{J} (\bar{t}) / μ^{J - 1}

⁠, in order to uniformly minimize the relative error for Eq. throughout the

\bar{t}

interval (the interested reader can resort to Ref. 38 for an example of an estimation for

J

⁠). Assuming now truncation at term

S_{j} (\bar{t}) / μ^{J - 1}

(J ⩽ J)

and since

\sum S_{j} (\bar{t}) / μ^{j - 1}

appears in the exponent of Eq. , it should hold for all

\bar{t}

that

| S_{j + 1} (\bar{t}) | / μ^{J} ≪ 1

in order for the WKB formula to be a good approximation for

a^{μ} (\bar{t})

⁠. In this case, the relative error between

a^{μ} (\bar{t})

and its associated WKB approximation becomes

\frac{a^{μ} (\bar{t}) - e^{μ \sum_{j = 0}^{J} S_{j} (\bar{t}) / μ^{j}}}{a^{μ} (\bar{t})} \sim 1 - O [e^{- S_{j + 1} (\bar{t}) / μ^{J}}] \sim O [S_{j + 1} (\bar{t}) / μ^{J}] .

(A23)

In summary, the expression of Eq. becomes useful if both of the following conditions are met:

| S_{j + 1} (\bar{t}) | < < μ | S_{j} (\bar{t}) |, | S_{j + 1} (\bar{t}) | < < μ^{J} .

(A24)

At this point, let us substitute the WKB ansatz of Eq. to Eq. and emphasize on the two leading-order terms of the WKB series,

\begin{align} \frac{1}{T^{2}} \{μ^{2} {[\frac{d S_{0} (\bar{t})}{d \bar{t}}]}^{2} + μ [\frac{d^{2} S_{0} (\bar{t})}{d {\bar{t}}^{2}} + 2 \frac{d S_{0} (\bar{t})}{d \bar{t}} \frac{d S_{1} (\bar{t})}{d \bar{t}}] + O (μ^{0})\} - \frac{1}{T} t r [M (\bar{t})] [μ \frac{d S_{0} (\bar{t})}{d \bar{t}} + \frac{d S_{1} (\bar{t})}{d \bar{t}} + O (μ^{- 1})] \\ - \frac{1}{T^{2}} / \frac{d m^{12} (\bar{t})}{d \bar{t}} m^{12} (\bar{t}) [μ \frac{d S_{0} (\bar{t})}{d \bar{t}} + O (μ^{0})] + \det [M (\bar{t})] + \frac{1}{T} m^{11} (\bar{t}) / \frac{d m^{12} (\bar{t})}{d \bar{t}} m^{12} (\bar{t}) - \frac{1}{T} \frac{d m^{11} (\bar{t})}{d \bar{t}} = 0 . \end{align}

(A25)

Among the T⁻² and T⁻¹ terms, the largest ones are

μ^{2} {[d S_{0} (\bar{t}) / d \bar{t}]}^{2}

and

μ tr [M (\bar{t})] d S_{0} (\bar{t}) / d \bar{t}

⁠, respectively (the reader should keep in mind that μ is taken to be a very large parameter). By dominant balance, the previous factors should have the same order of magnitude as

\det [M (\bar{t})]

(i.e., the coefficient of the T⁰ term), which in turn implies that O(μ²/T²) = O(μ/T) = O(1) and thus μ needs to be proportional to T (for simplicity, we choose μ = T). Of course, we have assumed in the above that

tr [M (\bar{t})] \neq 0

and

\det [M (\bar{t})] \neq 0

⁠, but an analogous rationale applies in these instances. {For the traceless and symmetric Hamiltonian model of Sec. II B where in general

\det [M (\bar{t})] \neq 0

⁠, it can be seen that O(μ²/T²) = O(1) and thus μ can be again taken to be equal to T.} After substituting μ = T in Eq. and taking the overall coefficients of T⁰ and T⁻¹ to be zero, the following formulas are acquired:

{[\frac{d S_{0} (\bar{t})}{d \bar{t}}]}^{2} - t r [M (\bar{t})] \frac{d S_{0} (\bar{t})}{d \bar{t}} + \det [M (\bar{t})] = 0,

(A26a)

\frac{d S_{1} (\bar{t})}{d \bar{t}} = - \frac{\frac{d^{2} S_{0} (\bar{t})}{d {\bar{t}}^{2}} - \frac{d m^{12} (\bar{t}) / d \bar{t}}{m^{12} (\bar{t})} [\frac{d S_{0} (\bar{t})}{d \bar{t}} - m^{11} (\bar{t})] - \frac{d m^{11} (\bar{t})}{d \bar{t}}}{2 \frac{d S_{0} (\bar{t})}{d \bar{t}} - t r [M (\bar{t})]} .

(A26b)

The first of the aforementioned relations represents the characteristic equation of matrix

M (\bar{t})

and thus accepts as solutions the eigenvalues

λ^{\pm} (\bar{t})

⁠,

\frac{d S_{0}^{\pm} (\bar{t})}{d \bar{t}} = λ^{\pm} (\bar{t}) \Leftrightarrow S_{0}^{\pm} (\bar{t}) = \int_{\bar{t}} λ^{\pm} ({\bar{t}}^{'}) d {\bar{t}}^{'} .

(A27)

Reexpressing the previous integral in terms of the temporal coordinate t, we obtain

S_{0}^{\pm} (t) = \frac{1}{T} \int_{t} λ^{\pm} (t^{'}) d t^{'} = \frac{1}{T} i ϕ^{\pm} (t) \Leftrightarrow T S_{0}^{\pm} (t) = i ϕ^{\pm} (t),

(A28)

with ϕ^±(t) standing for the standard dynamical phase [see also Appendix B 1, where ϕ^±(t) is attained via the mode population equations]. Given now that

d S_{0}^{\pm} (\bar{t}) / d \bar{t} = λ^{\pm} (\bar{t})

⁠, it can be seen for the denominator of Eq. that

\begin{align} 2 λ^{\pm} (\bar{t}) - t r [M (\bar{t})] & = λ^{\pm} (\bar{t}) - m^{11} (\bar{t}) + λ^{\pm} (\bar{t}) - m^{22} (\bar{t}) \\ = m^{12} (\bar{t}) r^{\pm} (\bar{t}) + \frac{m^{21} (\bar{t})}{r^{\pm} (\bar{t})} \\ = \frac{m^{12} (\bar{t}) {[r^{\pm} (\bar{t})]}^{2} + m^{21} (\bar{t})}{r^{\pm} (\bar{t})}, \end{align}

(A29)

where we have employed Eq. . [In the preceding set of equalities and in what follows,

r^{\pm} (\bar{t})

will signify

r_{R}^{\pm} (\bar{t})

as already indicated in Eq. .] Meanwhile, in the numerator of Eq. , we can substitute

λ^{\pm} (\bar{t}) = m^{12} (\bar{t}) r^{\pm} (\bar{t}) + m^{11} (\bar{t})

given Eq. and, subsequently, attain that

\begin{align} \frac{d S_{1}^{\pm} (\bar{t})}{d \bar{t}} & = - \frac{r^{\pm} (\bar{t}) d r^{\pm} (\bar{t}) / d \bar{t}}{{[r^{\pm} (\bar{t})]}^{2} + m^{21} (\bar{t}) / m^{12} (\bar{t})} \\ = - \frac{1}{2} \frac{d {[r^{\pm} (\bar{t})]}^{2} / d \bar{t}}{{[r^{\pm} (\bar{t})]}^{2} + m^{21} (\bar{t}) / m^{12} (\bar{t})} . \end{align}

(A30)

After integrating the latter equation, it comes that

\begin{align} S_{1}^{\pm} (\bar{t}) & = - \int_{\bar{t}} \frac{1}{2} \frac{d {[r^{\pm} ({\bar{t}}^{'})]}^{2} / d {\bar{t}}^{'}}{{[r^{\pm} ({\bar{t}}^{'})]}^{2} + m^{21} ({\bar{t}}^{'}) / m^{12} ({\bar{t}}^{'})} d {\bar{t}}^{'} \\ = - \int_{\bar{t}} \frac{1}{2} \frac{d {[r^{\pm} ({\bar{t}}^{'})]}^{2}}{{[r^{\pm} ({\bar{t}}^{'})]}^{2} + m^{21} ({\bar{t}}^{'}) / m^{12} ({\bar{t}}^{'})} . \end{align}

(A31)

The quantity in the integral represents the inner product

⟨V_{L}^{\pm} ({\bar{t}}^{'})| d V_{R}^{\pm} ({\bar{t}}^{'})⟩

based on Eqs. , , and . In other words, the following applies:

\begin{align} S_{1}^{\pm} (\bar{t}) = - \int_{\bar{t}} ⟨V_{L}^{\pm} ({\bar{t}}^{'})| d V_{R}^{\pm} ({\bar{t}}^{'})⟩ & = - \int_{\bar{t}} \frac{⟨v_{L}^{\pm} ({\bar{t}}^{'})| d v_{R}^{\pm} ({\bar{t}}^{'})⟩}{⟨v_{L}^{\pm} ({\bar{t}}^{'})| v_{R}^{\pm} ({\bar{t}}^{'}⟩} \\ = i γ_{G}^{\pm} (\bar{t}), \end{align}

(A32)

where

γ_{G}^{\pm} (\bar{t})

represents the geometric phase factor. [The form of

γ_{G}^{\pm} (\bar{t})

is also retrieved in Appendix B 1 in terms of the mode population equations, while the effect of the chosen eigenvector normalization is being extensively discussed in Appendix B 2.]

Putting it all together [Eqs. , , and ],

a^{T} (\bar{t})

or simply a(t) can be written as a linear superposition of the (±) solutions/modes,

a (t) \sim c_{a}^{+} e^{i γ_{G}^{+} (t)} e^{i ϕ^{+} (t)} + c_{a}^{-} e^{i γ_{G}^{-} (t)} e^{i ϕ^{-} (t)},

(A33)

where

c_{a}^{\pm}

are constants to be specified from initial conditions. The above formula is in agreement with Eq. (see Sec. III) after setting the integration limits from 0 to t and substituting

c_{a}^{\pm} (t) = c_{a}^{\pm} e^{i γ_{G}^{\pm} (t)}

⁠. {Had we considered higher-order WKB terms, then

c_{a}^{\pm}

would have been a function of time [i.e.,

c_{a}^{\pm} (t)

] to accommodate for such additional corrections and then

c_{a}^{\pm} (t) = c_{a}^{\pm} (t) e^{i γ_{G}^{\pm} (t)}

⁠.} According to Eq. (here, μ = T and J = 1), the following constraints must be met:

| S_{1}^{\pm} (\bar{t}) | < < T | S_{0}^{\pm} (\bar{t}) |, | S_{2}^{\pm} (\bar{t}) | < < T,

(A34)

for Eq. to be applicable, which then will differ from the exact solution by terms of order T⁻¹ given Eq. . Clearly, for sufficiently large T, both of the preceding conditions can be satisfied, unless

S_{1}^{\pm} (\bar{t}), S_{2}^{\pm} (\bar{t})

are unbounded [

S_{0}^{\pm} (\bar{t})

is bounded and this can be inferred from Eq. and the boundedness of the eigenvalues

λ^{\pm} (\bar{t})

]. Regarding function

S_{1}^{\pm} (\bar{t})

⁠, this cannot be predetermined solely based on Eq. or . The fact that the integrand quantity exhibits a singularity/pole when the Hamiltonian system operates at the spectrum degeneracy point [see Eq. in correlation with Eq. , or alternatively Eqs. , , and –, along with Eqs. and ] does not necessarily preclude that the result of the integration will display an analogous anomaly. (Consider, for instance, the integral of the inverse square root function

\int 1 / \sqrt{z} d z = 2 \sqrt{z} + c

⁠, with c being an additive constant: while the integrand becomes infinite at zero, the same is not true for the resulting integral expression.) It should be highlighted that for cyclic Hermitian settings, the aforementioned singularity essentially signifies the divergent behavior of the Berry connection at the DPs, which in turn has been associated with the emergence of topological phase transitions in periodic lattice systems^39,40 [i.e., in arrangements exhibiting a band structure, and where the reciprocal (or momentum) space plays the role of the parameter space]. Analogous observations have been reported in the presence of EPs for non-Hermitian settings, where an even wider range of nontrivial topological phases has been shown to be supported.^41,42 {Note that the Berry connection is, in general, complex-valued and will be given by

A^{\pm} (R) = i ⟨V_{L}^{\pm} (R)| \nabla_{R} |V_{R}^{\pm} (R)⟩

[this yields the Berry phase as

γ_{B}^{\pm} = \oint_{C^{R}} A^{\pm} (R) d R

⁠, given Eq. and identity

|d V_{R}^{\pm} (R)⟩ = \nabla_{R} |V_{R}^{\pm} (R)⟩ d R

with

\nabla_{R}

representing the nabla operator over

R

space], while under Hermitian conditions

(|V_{L}^{\pm} (R)⟩ = |V_{R}^{\pm} (R)⟩ = |V^{\pm} (R)⟩)

⁠, it comes that

A^{\pm} (R) = i ⟨V^{\pm} (R)| \nabla_{R} |V^{\pm} (R)⟩

and this can be shown to be real starting from the normalization condition

⟨V^{\pm} (R)| V^{\pm} (R)⟩ = 1

(see Appendix B 1 for more details). In the previous, we have employed the 1–1 correspondence between the scaled time variable

\bar{t}

and the trajectory

C^{R}

followed in parameter space

R

⁠, as described toward the end of the present subappendix.}

It is noteworthy to mention that all higher-order

S_{j}^{\pm} (\bar{t})

terms exhibit poles at spectrum degeneracies. To see this, we have considered in Eq. , after of course setting μ = T, the coefficients associated with factors T^−j (j ⩾ 2). The following difference–differential equation was then attained:

\begin{align} \frac{d^{2} S_{j - 1}^{\pm} (\bar{t})}{d {\bar{t}}^{2}} + \sum_{k = 0}^{j} \frac{d S_{k}^{\pm} (\bar{t})}{d \bar{t}} \frac{d S_{j - k}^{\pm} (\bar{t})}{d \bar{t}} - \frac{d m^{12} (\bar{t}) / d \bar{t}}{m^{12} (\bar{t})} \frac{d S_{j - 1}^{\pm} (\bar{t})}{d \bar{t}} \\ - tr [M (\bar{t})] \frac{d S_{j}^{\pm} (\bar{t})}{d \bar{t}} = 0, \end{align}

(A35)

which can be solved in terms of

d S_{j}^{\pm} (\bar{t}) / d \bar{t}

as

\frac{d S_{j}^{\pm} (\bar{t})}{d \bar{t}} = - \frac{\frac{d^{2} S_{j - 1}^{\pm} (\bar{t})}{d {\bar{t}}^{2}} + \sum_{k = 1}^{j - 1} \frac{d S_{k}^{\pm} (\bar{t})}{d \bar{t}} \frac{d S_{j - k}^{\pm} (\bar{t})}{d \bar{t}} - \frac{d m^{12} (\bar{t}) / d \bar{t}}{m^{12} (\bar{t})} \frac{d S_{j - 1}^{\pm} (\bar{t})}{d \bar{t}}}{2 \frac{d S_{0}^{\pm} (\bar{t})}{d \bar{t}} - t r [M (\bar{t})]} .

(A36)

The above expression reveals the divergence of the derivatives of

S_{j}^{\pm} (\bar{t})

for j ≥ 2 when

λ^{\pm} (\bar{t})

coalesce [in a similar fashion to

d S_{1}^{\pm} (\bar{t})

in Eq. ], and this can be seen from Eq. [λ_EP = tr(M_EP)/2] in conjunction with Eq.

[d S_{0}^{\pm} (\bar{t}) / d \bar{t} = λ^{\pm} (\bar{t})]

⁠. {We should emphasize here that the superscript ± in the higher-order

S_{j}^{\pm} (\bar{t})

terms stems from the presence of

S_{0}^{\pm} (\bar{t})

in their defining equation [i.e., Eq. ], and we know already that

S_{0}^{\pm} (\bar{t})

accepts two separate solutions associated with the eigenstates of the two-level Hamiltonian arrangement. In essence, each eigenmode is described by its own set of coefficients

S_{j}^{\pm} (\bar{t})

⁠, and had we desired a higher-order WKB approximation, a(t) could have been written as a linear superposition of factors

e^{\sum S_{j}^{\pm} (\bar{t}) / T^{j - 1}}

bearing in mind the constraints of Eq. along with equality μ = T.}

Having provided the conditions [Eq. ] for the WKB result of Eq. to be a satisfactory approximation, we shall now compute the factor

S_{1}^{\pm} (\bar{t})

in the special case where

m^{12} (\bar{t}) = m^{21} (\bar{t})

[i.e., matrix

M (\bar{t})

—or equivalently the Hamiltonian

H (\bar{t})

—is symmetric]. The integration in Eq. can be then directly performed and will yield

S_{1}^{\pm} (\bar{t}) = - \frac{1}{2} \ln ({[r^{\pm} (\bar{t})]}^{2} + 1) + c,

(A37)

where c is a constant that can be simply absorbed in factors

c_{a}^{\pm}

of Eq. . If it further applies that

M (\bar{t})

[and similarly

H (\bar{t})

] is traceless, then the non-Hermitian model of Sec. II B is attained. From Eq. , it can be subsequently deduced that

2 λ^{\pm} (\bar{t}) = m^{12} (\bar{t}) \frac{{[r^{\pm} (\bar{t})]}^{2} + 1}{r^{\pm} (\bar{t})},

(A38)

which leads to an alternative expression to Eq. ,

S_{1}^{\pm} (\bar{t}) = - \frac{1}{2} \ln (\frac{2 λ^{\pm} (\bar{t}) r^{\pm} (\bar{t})}{m^{12} (\bar{t})}) + c .

(A39)

It can be seen from either Eq. or that the conditions depicted in Eq. cannot be satisfied at the EP degeneracies, as

S_{1}^{\pm} (\bar{t})

becomes infinite there [

r_{EP}^{2} = - 1

from Eq. after substituting

m_{EP}^{12} = m_{EP}^{21}

⁠, and additionally λ_EP = 0 for traceless systems given Eq. ]. In other words, at least for the special case of symmetric

M (\bar{t})

or

H (\bar{t})

matrices, the WKB approach is not applicable at the location of EPs (the traceless property does not affect this conclusion). In Appendix B 3, where the mode population equations are analytically investigated under non-Hermitian conditions, it will be shown how the latter shortcoming can be resolved by studying the limiting behavior of the quantum arrangement at the vicinity of such degeneracies.

To summarize, we shall provide the following concluding notes:

In the present appendix, we have shown how via the WKB method an approximate expression can be attained for the state vector component a(t) in the adiabatic limit T → ∞ [analogous formulas can also be derived for component b(t)]. A direct correlation between the two leading-order terms of the WKB series $\sum S_{j}^{\pm} (\bar{t}) / T^{j - 1}$ with the dynamical [ϕ^±(t)] and geometric $[γ_{G}^{\pm} (\bar{t})]$ phases has been demonstrated [see Eqs. (A28) and (A32)]. The conditions for approximating a(t) as the product of the aforementioned phase factors [see Eq. (A33)] have been outlined [see Eq. (A34)], while in the same time, a difference–differential relation has been acquired that can yield all higher-order terms/phases $S_{j}^{\pm} (\bar{t}) / T^{j - 1}$ with j ≥ 2 [see Eq. (A35) along with the general constraints for the validity of the WKB method in Eq. (A24)]. It should be highlighted that exactly the same results can be attained had we originally emphasized on the differential equation that the transformed component $\tilde{a} (t)$ satisfies [see Eq. (12) of Sec. II A], after of course reinstating a(t) in the retrieved WKB formulas [see Eq. (10) of Sec. II A].
A few notes must also be made regarding the terminology. In Appendix A 1, the system matrix M (and Hamiltonian H) and its eigenvalues and eigenvectors (and the associated eigenvector component ratios) were written as functions of time t. As explained though in Appendix A 1 (and at the beginning of Sec. II A), the above quantities depend directly on the scaled time coordinate $\bar{t}$ ⁠, and that is why, here, they were written as a function of $\bar{t}$ (the same applies also for terms $S_{j}^{\pm}$ ⁠, including the geometric phase $γ_{G}^{\pm}$ ⁠). This comes in contrast with the state vector component a (and b), which depends on both $\bar{t}$ and T as becomes readily evident from Eqs. (A21) and (A22) (the latter in essence implies that if t and T are multiplied by the same factor, a will vary despite that $\bar{t}$ remains unchanged). Consequently, we write $a^{T} (\bar{t})$ [or simply a(t)] and not $a (\bar{t})$ ⁠, to correctly reflect the underlying functional dependencies. In a similar fashion, the dynamical phase ϕ^± is a function of both $\bar{t}$ and T [see Eqs. (A27) and (A28)] and thus can be written equivalently as ϕ^±(t) or $ϕ_{\pm}^{T} (\bar{t})$ [in the analysis that follows, we might also use the auxiliary variable $Φ^{\pm} (\bar{t}) = ϕ^{\pm} (t) / T$ ⁠, which will depend solely on $\bar{t}$ ]. To summarize, an arbitrary quantity q that is a function only of $\bar{t}$ will be written alternatively as q(t) or $q (\bar{t})$ [see for instance $S_{0}^{\pm}$ or λ^± in Eqs. (A27) and (A28)], while if q is a function of both $\bar{t}$ and T, then this will be expressed as q(t) or $q^{T} (\bar{t})$ [see the component a in Eqs. (A21), (A22), and (A33)].
A fundamental difference between the dynamical and geometric phases can be underlined at this point: the former depends on how fast or slow the Hamiltonian varies with time (i.e., the magnitude of T), while the latter depends solely on the geometry of the path $C^{R}$ followed within parameter space $R$ (thus, the term “geometric phase”). Essentially, each point along $C^{R}$ can be described by a specific value of $\bar{t}$ ⁠, which is independent from the choice of period T (this does not apply for t, as if for example T is increased, then t will signify points along $C^{R}$ that are closer to the starting point of evolution). This stems from the fact that the various Hamiltonian parameters are functions only of $\bar{t}$ in our adiabatic analysis, and thus, $C^{R}$ can be fully parameterized in terms of $\bar{t}$ (and not of t). Consequently, all points of $C^{R}$ will be characterized by a unique $γ_{G}^{\pm} (\bar{t})$ (the 1–1 correspondence is guaranteed via variable $\bar{t}$ ⁠), but not from a unique ϕ^±(t) (given that the latter depends not only on $\bar{t}$ but also on T). Moreover, it should be noted that given the 1–1 mapping between $\bar{t}$ and $C^{R}$ ⁠, the integrals in Eqs. (A27), (A31), and (A32) can be equivalently expressed as line integrals along $C^{R}$ ⁠.

3. Asymptotic behavior of Airy functions for arguments of large absolute value

By this point, the efficacy of the WKB approach in analyzing the slowly varying dynamics of generalized Hamiltonian settings within high orders of accuracy has been clearly illustrated. In the present and following sections, we shall apply such a methodology to attain asymptotic expressions for two distinct classes of functions, namely the Airy ( Appendix A 3) and the modified Bessel ( Appendix A 4) functions.

Let us first state the Airy equation, which reads as

\frac{d^{2} w (z)}{d z^{2}} - z w (z) = 0,

(A40)

and is defined over the complex plane. To make use of the WKB method, Eq. is reexpressed as follows:

\frac{1}{M^{2}} \frac{d^{2} W^{μ} (z)}{d z^{2}} - z W^{μ} (z) = 0,

(A41)

where

M

is assumed to be a very large positive real quantity (this is a standard trick in asymptotic analysis, and for more examples, one can resort to Ref. 37). After substituting the WKB ansatz

(μ \in R^{+}, μ \to \infty)

W^{μ} (z) \sim e^{μ \sum_{j = 0}^{\infty} S_{j} (z) / μ^{j}}

(A42)

in Eq. , it comes that

\begin{align} \frac{1}{M^{2}} [μ^{2} {(\frac{d S_{0}}{d z})}^{2} + μ (\frac{d^{2} S_{0}}{d z^{2}} + 2 \frac{d S_{0}}{d z} \frac{d S_{1}}{d z}) + \frac{d^{2} S_{1}}{d z^{2}} \\ + 2 \frac{d S_{0}}{d z} \frac{d S_{2}}{d z} + {(\frac{d S_{1}}{d z})}^{2} + O (μ^{- 1})] = z, \end{align}

(A43)

where we have omitted the dependencies from z to simplify our notation. The largest term on the left-hand side is

μ^{2} {(d S_{0} / d z)}^{2} / M^{2}

⁠, which by dominant balance should be of an equal order of magnitude as z appearing on the right-hand side. In other words, μ should be proportional to

M

and we can simply choose

μ = M

⁠. By taking now the coefficients of

M^{0}

⁠,

M^{- 1}

⁠, and

M^{- 2}

in Eq. to be equal to zero (after of course setting

μ = M

⁠), the following relations are, respectively, acquired:

\frac{d S_{0}}{d z} = \pm \sqrt{z} \Leftrightarrow S_{0} = \pm \frac{2}{3} z^{3 / 2} + c,

(A44a)

\begin{gathered} \frac{d S_{1}}{d z} = - \frac{1}{2} \frac{d^{2} S_{0} / d z^{2}}{d S_{0} / d z} \Leftrightarrow \\ S_{1} = - \frac{1}{2} \ln (\frac{d S_{0}}{d z}) + c = - \frac{1}{4} \ln (z) + c^{'}, \end{gathered}

(A44b)

\begin{gathered} \frac{d^{2} S_{1}}{d z^{2}} + 2 \frac{d S_{0}}{d z} \frac{d S_{2}}{d z} + {(\frac{d S_{1}}{d z})}^{2} = 0 \Leftrightarrow \\ \frac{d S_{2}}{d z} = \mp \frac{5}{32} z^{- 5 / 2} \Leftrightarrow S_{2} = \pm \frac{5}{48} z^{- 3 / 2} + c . \end{gathered}

(A44c)

[There is a ± sign emerging from dS₀/dz in Eq. , which has been absorbed in constant c′.] If we further assume truncation of series

\sum S_{j} (z) / M^{j - 1}

at the second term (i.e., at j = 1), then for the WKB method to yield a valid approximation, it must apply that

| S_{1} | < < M | S_{0} |, | S_{2} | < < M,

(A45)

according to the discussion in Appendix A 2 [see Eq. or ]. The above constraints need not be examined for any value of

M

⁠, but only for

M = 1

⁠, which describes Eq. . Given also Eq. , we attain

| \ln (z) | < < | z |^{3 / 2}, | z |^{- 3 / 2} < < 1,

(A46)

which both apply for sufficiently large magnitudes of z. Based then on Eqs. and , the leading-order asymptotic behavior of the Airy functions when |z| → ∞ becomes

w (z) \sim c_{1} \frac{e^{2 z^{3 / 2} / 3}}{\sqrt[4]{z}} + c_{2} \frac{e^{- 2 z^{3 / 2} / 3}}{\sqrt[4]{z}},

(A47)

with c_1,2 constants to be determined from boundary conditions, connection formulas, or other restrictions/information particular to the problem at hand. For instance, if we desire to attain the exact asymptotic expression for Ai(z) (Airy function of the first kind), we have to first consider that when

z \in R^{+}

⁠, Ai(z → ∞) = 0. Consequently, c₁ = 0 at least when z is a positive real, so as Eq. does not blow up. The exact value for c₂ can come from the integral representation of Ai(z) for

z \in R

⁠,³¹^,

A i (z) = \frac{1}{2 π} \int_{- \infty}^{\infty} e^{i x^{3} / 3 + i z x} d x = \frac{1}{π} \int_{0}^{\infty} \cos (\frac{x^{3}}{3} + z x) d x,

(A48)

where after applying the method of stationary phases (known also as steepest descent in complex analysis) for z → ∞, it comes that⁴³

A i (z) \sim \frac{1}{2 \sqrt{π}} \frac{e^{- 2 z^{3 / 2} / 3}}{\sqrt[4]{z}},

(A49)

and thus

c_{2} = {(2 \sqrt{π})}^{- 1}

⁠. [We should highlight here, based also on Eq. , that Ai(z) attains real values for real z, and this should come as no surprise since Eq. supports real solutions under such circumstances.] Alternatively, this can be deduced from the following relations:

A i (z) = \frac{1}{π} \sqrt{\frac{z}{3}} K_{1 / 3} (\frac{2}{3} z^{3 / 2}),

(A50a)

K_{1 / 3} (z) \sim \sqrt{\frac{π}{2 z}} e^{- z},

(A50b)

which both apply when

z \in R^{+}

and the latter when z → ∞³¹ [K_1/3(·) refers to the modified Bessel function of the second kind of order 1/3].

Of interest is now to notice that when z < 0, the solutions to Eq. will exhibit an oscillatory behavior. {For z > 0, Eq. accepts exponentially decaying or growing solutions, which is readily reflected in the dynamics of the Airy functions of the first [see Eq. ] and second (see Ref. 31) kind.} Indeed, Ai(z) oscillates for negative values of z and keeping in mind that

A i (z) \in R

when

z \in R

[Eq. ], it comes from Eq. that if z approaches the negative real axis from the second (third) quadrant, i.e., arg(z) → π [arg(z) → −π], it must be true that

c_{1} = i c_{2} = i {(2 \sqrt{π})}^{- 1}

[c_{1} = - i c_{2} = - i {(2 \sqrt{π})}^{- 1}]

given that the principal values are employed for arg(z) [−π < arg(z) ⩽ π]. As a result, the following formula applies for large negative z = ρe^iπ = −ρ or positive ρ:

A i (- ρ) \sim \frac{1}{\sqrt{π}} ρ^{- 1 / 4} \sin (\frac{2}{3} ρ^{3 / 2} + \frac{π}{4}),

(A51)

and this can be independently verified by applying the method of stationary phases on Eq. for z → −∞⁴³ [or alternatively by expressing Ai(z), when z < 0, in terms of Bessel functions and subsequently exploit the very well-known asymptotic formulas of the latter³¹].

At this point, we can make an interesting observation: the coefficient c₁ changes from zero when z > 0, to a nonzero value when z < 0. So how can it be that Ai(z) is analytic everywhere, but it does not have a uniform asymptotic expansion throughout the complex plane? This famous paradox in complex analysis is known as the Stokes phenomenon and is largely attributed to the presence of multivalued functions in the asymptotic expressions of analytic functions. [Consider, for instance, $e^{2 z^{3 / 2} / 3} z^{- 1 / 4}$ and that we employ the principal branches for both (multivalued) functions z^3/2 and z^−1/4 (branch cut is then along axis z < 0). Taking |z| = ρ > 0, it comes that as arg(z) → π, then $e^{2 z^{3 / 2} / 3} z^{- 1 / 4} \to e^{- 2 i ρ^{3 / 2} / 3 - i π / 4} ρ^{- 1 / 4}$ ⁠, while as arg(z) → −π, then $e^{2 z^{3 / 2} / 3} z^{- 1 / 4} \to e^{2 i ρ^{3 / 2} / 3 + i π / 4} ρ^{- 1 / 4}$ ⁠. Hence, at the same point along the negative real axis (z = −ρ), $e^{2 z^{3 / 2} / 3} z^{- 1 / 4}$ attains different values and this dictates that c₁ should have different values below and above axis z < 0. Such a fact has been already pointed out when deriving the asymptotic formula of Eq. (A51).]

To resolve the aforementioned inconsistency, Stokes considered what happens as z moves along a circle around the origin of the complex plane. He realized that the exponential factors in Eq. (A47) assumed both large and small absolute values intermittently: while one term is superior, the other is inferior and vice versa (this is contingent upon whether the exponent has a positive or a negative real part). Furthermore, a “jump” in the value of a coefficient (c₁ or c₂) must not compromise the integrity of the asymptotic approximation as a whole—if a new exponential term is set to emerge due to the switching of its coefficient from zero to nonzero, it should occur when it is the least noticeable. This led to Stokes’ ingenious observation that a jump could only manifest at points where the switching coefficient in question corresponds to an inferior term and the superior term reaches its maximum.⁴⁴ (In a much later study,⁴⁵ Berry actually showed that such abrupt jumps do not occur discontinuously but rather in a smooth fashion, by truncating the dominant part of the asymptotic series expansion near its smallest element and subsequently studying the subdominant contributions via Dingle’s method of Borel sums.⁴⁶) These points are known to lie on the Stokes lines, which in Eq. (A47) occur when Im(z^3/2) = 0 ⇔ arg(z) = 0, ±2π/3. [Note that the locations where Re(z^3/2) = 0 ⇔ arg(z) = ±π/3, π define the anti-Stokes lines, and there $e^{2 z^{3 / 2} / 3}$ and $e^{- 2 z^{3 / 2} / 3}$ become comparable in magnitude.] Then, for arg(z) = 0, it is $e^{2 z^{3 / 2} / 3}$ that asymptotically dominates over $e^{- 2 z^{3 / 2} / 3}$ ⁠, while this is reversed when arg(z) = ±2π/3. Hence, if c₁ (c₂) exhibits a jump, this better takes place when arg(z) = ±2π/3 [arg(z) = 0]. The above rationale is of very general application to a class of functions frequently emerging in physical problems and exhibiting rapidly varying (oscillating or exponentially growing/decaying) components (examples include the Bessel, parabolic cylinder, and hypergeometric functions, or even the Fresnel and elliptic integrals). In this respect, if the leading behaviors in Eq. (A47) are given by e^P(z) and e^Q(z) (which are also assumed to display significantly different dynamics across the complex space), the (anti-)Stokes curves will be described by the equation Im[P(z) − Q(z)] = 0 {Re[P(z) − Q(z)] = 0}, and in their vicinity, one of the solutions will prevail (both solutions will be comparable).⁴³

Returning to the specific case of the Airy function of the first kind Ai(z), given that for

z \in R^{+}

⁠, it becomes exponentially decaying, we need not concern about the Stokes line at arg(z) = 0. [This comes in contrast to the Airy function of the second kind Bi(z), which is exponentially growing for positive z and undergoes a discontinuous change of c₂ at arg(z) = 0.] Consequently, we emphasize on the Stokes lines at arg(z) = ±2π/3, where c₁ changes from zero to a non-zero value, enabling thus factor

e^{2 z^{3 / 2} / 3}

to come into play in Eq. . From the previous discussion, it is expected that Eq. will remain applicable to arg(z) in the range −2π/3 to 2π/3. Taking advantage of the connection formula,

e^{i 2 π / 3} A i (z^{'} e^{i 2 π / 3}) + e^{- i 2 π / 3} A i (z^{'} e^{- i 2 π / 3}) + A i (z^{'}) = 0,

(A52)

and assuming that 0 < arg(z′) ⩽ π/3 [−π/3 < arg(z′) ⩽ 0], the asymptotic representation of Ai(z) can be retrieved in the domain 2π/3 < arg(z) ⩽ π [−π < arg(z) ⩽ −2π/3] after replacing Ai(z′), Ai(z′e^−i2π/3) [Ai(z′e^i2π/3)] by Eq. and subsequently setting z = z′e^i2π/3 (z = z′e^−i2π/3). It was then found that

c_{1} = i c_{2} = i {(2 \sqrt{π})}^{- 1}

[c_{1} = - i c_{2} = - i {(2 \sqrt{π})}^{- 1}]

⁠, which in turn confirms our asymptotic results for z < 0 as arg(z) → π [arg(z) → −π]. [For the latter, see Eq. and comments above.]

We are in the position now to state the leading-order asymptotics of the Airy function of the first kind throughout the complex plane,

A i (z) \sim \{\begin{matrix} - \frac{i e^{\frac{2 z^{3 / 2}}{3}}}{2 \sqrt{π} \sqrt[4]{z}} + \frac{e^{- \frac{2 z^{3 / 2}}{3}}}{2 \sqrt{π} \sqrt[4]{z}}, & \arg (z) ⩽ - \frac{2 π}{3}, \\ \frac{e^{- \frac{2 z^{3 / 2}}{3}}}{2 \sqrt{π} \sqrt[4]{z}}, & - \frac{2 π}{3} < \arg (z) ⩽ \frac{2 π}{3}, \\ \frac{i e^{\frac{2 z^{3 / 2}}{3}}}{2 \sqrt{π} \sqrt[4]{z}} + \frac{e^{- \frac{2 z^{3 / 2}}{3}}}{2 \sqrt{π} \sqrt[4]{z}}, & \arg (z) > \frac{2 π}{3}, \end{matrix}

(A53)

which agrees with the literature³¹ and moreover exhibits a relative error of O(z^−3/2) with respect to the exact form of Ai(z), as can be seen via Eqs. and [simply substitute

μ = M = 1

and J = 1 in Eq. and use the expression for S₂ from Eq. ]. It is noteworthy to mention that the leading approximations for dAi(z)/dz, Bi(z), dBi(z)/dz can be obtained from Eq. through direct differentiation [dAi(z)/dz] and by utilizing appropriate connection formulas [Bi(z), dBi(z)/dz]. [Directly differentiating Eq. , to attain the approximate form of dAi(z)/dz and subsequently of dBi(z)/dz for large z, has been well established by Theorem A in Ref. 47.] Alternatively, the methodology developed in the current section can be employed to acquire the asymptotic profiles of the preceding functions from scratch.

4. Asymptotic behavior of modified Bessel functions for orders of large absolute value

As an additional example, the approximate expressions for functions I_ν(νz), K_ν(νz) for large |ν| and finite z shall be derived (ν, z are in general complex). We will start with the differential equation that

\sqrt{z} I_{ν} (ν z), \sqrt{z} K_{ν} (ν z)

satisfy

\frac{1}{ν^{2}} \frac{d^{2} w^{μ} (z)}{d z^{2}} - (\frac{1 + z^{2}}{z^{2}} - \frac{1}{4 ν^{2} z^{2}}) w^{μ} (z) = 0 .

(A54)

Substituting the WKB formula

w^{μ} (z) \sim e^{μ \sum_{j = 0}^{\infty} S_{j} (z) / μ^{j}},

(A55)

it comes that

\begin{align} \frac{1}{v^{2}} [μ^{2} {(\frac{d S_{0}}{d z})}^{2} + μ (\frac{d^{2} S_{0}}{d z^{2}} + 2 \frac{d S_{0}}{d z} \frac{d S_{1}}{d z}) + \frac{d^{2} S_{1}}{d z^{2}} \\ + 2 \frac{d S_{0}}{d z} \frac{d S_{2}}{d z} + {(\frac{d S_{1}}{d z})}^{2} + \frac{1}{4 z^{2}} + O (μ^{- 1})] = \frac{1 + z^{2}}{z^{2}} . \end{align}

(A56)

As opposed to our former WKB calculations, here, we let μ to be complex and furthermore assume |μ| → ∞. Then, from dominant balance considerations, μ must have the same order of magnitude with v and we can simply set μ = v in the previous equation. After taking the coefficients of v⁰, v⁻¹, v⁻² to be zero, the expressions for the leading coefficients of the WKB series are computed as

\begin{align} \frac{d S_{0}}{d z} & = \pm \frac{\sqrt{1 + z^{2}}}{z} \Leftrightarrow S_{0} = \pm \sqrt{1 + z^{2}} \mp \sinh^{- 1} (\frac{1}{z}) + c \\ \Leftrightarrow S_{0} = \pm \sqrt{1 + z^{2}} \pm \ln (\frac{z}{1 + \sqrt{1 + z^{2}}}) + c = \pm γ + c, \end{align}

(A57a)

S_{1} = - \frac{1}{2} \ln (\frac{d S_{0}}{d z}) + c = - \frac{1}{2} \ln (\frac{\sqrt{1 + z^{2}}}{z}) + c^{'},

(A57b)

\begin{gathered} \frac{d S_{2}}{d z} = \pm \frac{4 z - z^{3}}{8} {(1 + z^{2})}^{- 5 / 2} \Leftrightarrow \\ S_{2} = \pm [\frac{1}{8} {(1 + z^{2})}^{- 1 / 2} - \frac{5}{24} {(1 + z^{2})}^{- 3 / 2}] + c . \end{gathered}

(A57c)

For the WKB analysis to produce an accurate estimate, the following conditions must be met:

| S_{1} | < < | ν S_{0} |, | S_{2} | < < | ν | .

(A58)

[The latter comes from Eq. or after replacing μ with |ν| and J with unity.] For sufficiently large |ν|, both of the aforementioned relations hold true, except when z = 0, ±i. These are identified as transition points [according to theory,³³ these correspond to the zeros of (1 + z²)/z (known also as turning points) and the singularities of (1 + z²)/z and 1/(4z²) in Eq. ] and signify the locations where terms S_j (j = 0, 1, 2) become unbounded. Apart from such inflection sites, the leading asymptotic form of the solution to Eq. can be described as

w^{ν} (z) \sim c_{1} e^{ν γ} \frac{\sqrt{z}}{\sqrt[4]{1 + z^{2}}} + c_{2} e^{- ν γ} \frac{\sqrt{z}}{\sqrt[4]{1 + z^{2}}},

(A59)

with a relative error O(ν⁻¹) [see Eq. and set μ = ν and J = 1]. The constants c_1,2 can be cleverly found by employing the familiar asymptotics of the modified Bessel functions,

\sqrt{z} I_{ν} (ν z) \sim \frac{1}{\sqrt{2 π ν}} e^{ν z}, \sqrt{z} K_{ν} (ν z) \sim \sqrt{\frac{π}{2 ν}} e^{- ν z},

(A60)

where z is taken to be large and positive, ν is considered to be fixed, and |arg(νz)| < π/2 ⇔ |arg(ν)| < π/2.³² Capitalizing on the uniform property of Eq. , it can be deduced from Eq. that

c_{1} = 1 / \sqrt{2 π ν}, c_{2} = \sqrt{π / (2 ν)}

[it should be reminded that Eq. admits the solutions

\sqrt{z} I_{ν} (ν z)

and

\sqrt{z} K_{ν} (ν z)

]. In the most general scenario then (z finite or infinite), it will be true as |ν| → ∞ that

I_{ν} (ν z) \sim \frac{1}{\sqrt{2 π ν}} \frac{e^{ν γ}}{\sqrt[4]{1 + z^{2}}}, K_{ν} (ν z) \sim \sqrt{\frac{π}{2 ν}} \frac{e^{- ν γ}}{\sqrt[4]{1 + z^{2}}},

(A61)

with |arg(ν)| < π/2 and |arg(νγ)| < π/2.³³ Clearly, Eqs. and are consistent, meaning that the latter is the limit of the former as z goes to plus infinity [from Eq. , γ → z and

\sqrt{z} / \sqrt[4]{1 + z^{2}} \to 1

⁠, as z grows large and positive]. If now ν or νγ lies in the left side of the complex plane, we need to resort to the subsequent connection,³²^,

I_{- ν^{'}} (ν^{'} z^{'}) = I_{ν^{'}} (ν^{'} z^{'}) + 2 π^{- 1} \sin (π ν^{'}) K_{ν^{'}} (z^{'}),

(A62a)

K_{- ν^{'}} (ν^{'} z^{'}) = K_{ν^{'}} (ν^{'} z^{'}),

(A62b)

and analytic continuation formulas,³²^,

I_{ν^{'}} (ν^{'} z^{'} e^{i k π}) = e^{i k π ν^{'}} I_{ν^{'}} (ν^{'} z^{'}),

(A63a)

\begin{align} K_{ν^{'}} (ν^{'} z^{'} e^{i k π}) & = - i π \sin (k π ν^{'}) \csc (π ν^{'}) I_{ν^{'}} (ν^{'} z^{'}) \\ + e^{- i k π ν^{'}} K_{ν^{'}} (ν^{'} z^{'}), \end{align}

(A63b)

with

k \in Z

(k = ±1 when using principal branches) and csc(·) representing the cosecant or inverse sine function [csc(·) = 1/sin(·)]. Equations and reveal that the expansions of I_ν(νz), K_ν(νz) are in general mixtures of positive and negative exponential-type solutions, while the leading-order forms of their derivatives can be attained simply by differentiating the respective asymptotic expressions.⁴⁷ Note that the specific case arg(ν) = ±π/2 requires special attention and has been covered extensively in Refs. 48 and 49.

At this point, we should highlight that great care must be given in all the aforementioned computations, owing to the presence of multivalued functions in the various formulas. The selection of a specific branch causes branch cuts to appear in complex space, which subsequently introduces complications even in the simplest of the calculations as laboriously depicted in Sec. II B. Yet, the effect of branch cuts has not been demonstrated so far in complex integration. In Appendixes A 2– A 4, we have avoided integrating along specific paths (indefinite integrals have been used in all instances), so as not to consider scenarios where such paths intersect with branch lines. [Had this occurred, the integral in question should have been evaluated separately before and after these discontinuity curves, as we shall subsequently see, and this in turn can give rise to additional multiplication constants in the WKB expressions of Eqs. (A33), (A47), and (A59), thus modifying the already existing constants $c_{a}^{\pm}, c_{1, 2}$ —the previously described mathematics are in essence the culprit behind the manifestation of the Stokes phenomenon (see also the detailed description of Appendix A 3)]. This simplifies the underlying analysis and enables us to better illustrate the utility of the WKB method through the examined asymptotic studies. In what follows, it will be clarified via certain examples how the path of integration affects the values of the associated complex integrals, particularly when branch lines are crossed or when the integrand displays singularities within the domain of integration.

5. A few notes on the evaluation of complex integrals in the presence of branch cuts and singularity points

Throughout the current investigation, we have employed for our calculations the principal values for the arguments ϑ of complex numbers, i.e., −π < ϑ ⩽ π (see the beginning of Sec. II B 1). Of course, the argument of a complex quantity is simply one of a garden variety of multivalued complex functions (examples are the logarithmic, the square root, the nth root with n being a positive integer, and the inverse trigonometric). In order to be mathematically consistent, we have chosen all such functions to be represented by their principal branches (for instance, $\sqrt{z} = \sqrt{ρ} e^{i ϑ / 2}$ ⁠, with z = ρe^iϑ, ρ > 0, −π < ϑ ⩽ π), which in turn implies the emergence of branch cuts in the complex plane where discontinuities arise (e.g., for the principal square/nth root and logarithm, the branch cut extends from the origin along the negative real axis to negative infinity). In general, branch cuts designate the locations that the different sheets/branches of a multivalued function meet and thus are typically employed in complex analysis to delimit the function’s domain in regions where it becomes single-valued and well-behaved (here, such a region corresponds to the principal branch).

Now, the question comes as to what happens if a line integral crosses such a branch cut? So far, in the evaluation of the various complex integrals, the latter scenario was not taken into account, as our emphasis was on demonstrating the effectiveness of the WKB method in analyzing the asymptotic behavior of generalized Hamiltonian settings (and subsequently of particular function categories), without delving into intricate mathematical details. In what follows, we shall show via specific examples how to address such special situations, while the effect of poles/singularities will also be separately discussed.

Let us first consider the integration of function 1/z, with z being a complex variable, along a non-self-intersecting, continuous, and differentiable (known also as simple and smooth) curve C with initial and final points at

z_{ι} = e^{i ϑ_{ι}} = e^{i 2 π / 3}

and

z_{f} = e^{i ϑ_{f}} = e^{- i 2 π / 3}

⁠, respectively. Moreover, we assume that C intersects the negative real axis (ϑ = π) once, at point

z_{b} = e^{i θ_{b}} = e^{i π} = - 1

⁠. Now, we have that

\begin{align} \int_{C} z^{- 1} d z & = \int_{C^{+}} z^{- 1} d z + \int_{C^{-}} z^{- 1} d z = \ln (z_{f}) - \ln (z_{b}^{+}) \\ + \ln (z_{b}^{-}) - \ln (z_{ι}) = i (ϑ_{f} - ϑ_{ι}) + i (ϑ_{b}^{-} - ϑ_{b}^{+}) = \frac{2 i π}{3}, \end{align}

(A64)

where curves C^± refer to the part of C just after z_b (C⁺:

z_{b}^{+} \to z_{f}

⁠, with

z_{b}^{+} = e^{i θ_{b}^{+}}

and

θ_{b}^{+}

signifying angles just larger than −π, i.e.,

θ_{b}^{+} = - π^{+}

⁠) and just before z_b (C⁻:

z_{ι} \to z_{b}^{-}

⁠, with

z_{b}^{-} = e^{i θ_{b}^{-}}

and

θ_{b}^{-}

signifying angles just smaller than π, i.e.,

θ_{b}^{-} = π^{-}

⁠). Note that the separate evaluation along C⁺ and C⁻ of the contour integral in Eq. was necessitated from the fact that the antiderivative ln(z) of 1/z is not defined along axis ϑ = π. To further verify our result, we can consider the integration of 1/z along a closed and positively oriented curve C′, which encircles the origin of the complex plane and intersects with ϑ = π again at z_b = −1. Then, given Eq. and allowing for the beginning (z_ι) and ending (z_f) points to be identical, it comes that

\oint_{C^{'}} z^{- 1} d z = i (ϑ_{b}^{-} - ϑ_{b}^{+}) = 2 i π .

(A65)

The above can be independently verified from the residue theorem of complex analysis. (Since 1/z is analytic everywhere except from an isolated point at the origin, the residue theorem is applicable.) Indeed, given that z = 0 is a simple pole of 1/z, it will apply that Res(1/z, 0) = lim_z→0 z · 1/z = 1, and then, the integral in Eq. should be equal to 2iπ Res(1/z, 0) or 2iπ (the winding number of C′ around z = 0 is equal to one, assuming a positive/counterclockwise orientation for C′).

Another nontrivial integration example is that of the square root function, which instead of an anomaly point exhibits a branch cut along the non-positive real axis. [This typically appears in the adiabatic analysis of two-mode Hamiltonian systems, owing to the presence of square roots in the formulas for the eigenvalues (and consequently of eigenvectors—see Appendix A 1), which are subsequently integrated to attain the various phase factors (see Appendix A 2)]. Following an analogous methodology as in Eq. , we can find that

\begin{align} \int_{C} \sqrt{z} d z & = \int_{C^{+}} \sqrt{z} d z + \int_{C^{-}} \sqrt{z} d z \\ = \frac{2}{3} [z_{f}^{3 / 2} - z_{ι}^{3 / 2} + {(z_{b}^{-})}^{3 / 2} - {(z_{b}^{+})}^{3 / 2}] = - \frac{4 i}{3} . \end{align}

(A66)

In the case of the closed and positively oriented path C′, it can be attained that

\oint_{C^{'}} \sqrt{z} d z = \frac{2}{3} [{(z_{b}^{-})}^{3 / 2} - {(z_{b}^{+})}^{3 / 2}] = - \frac{4 i}{3} .

(A67)

Of course, the residue theorem is no longer applicable since

\sqrt{z}

is nonanalytic over the whole axis ϑ = π instead of a finite number of anomaly points (axis ϑ = π is partially contained within the chosen path C′).

It is of interest to note that in the special case that C represents a segment of the unit circle and C′ the unit circle itself, all the aforementioned results could have been alternatively retrieved by employing the parameterization z = e^iϑ. More specifically, it would apply

\begin{align} \int_{C} z^{- 1} d z & = \int_{C^{+}} z^{- 1} d z + \int_{C^{-}} z^{- 1} d z \\ = i \int_{ϑ_{b}^{+}}^{ϑ_{f}} d ϑ + i \int_{ϑ_{ι}}^{ϑ_{b}^{-}} d ϑ \\ = i (ϑ_{f} - ϑ_{ι}) + i (ϑ_{b}^{-} - ϑ_{b}^{+}) = \frac{2 i π}{3}, \end{align}

(A68a)

\oint_{C^{'}} z^{- 1} d z = i (ϑ_{b}^{-} - ϑ_{b}^{+}) = 2 i π,

(A68b)

\begin{gathered} \int_{C} \sqrt{z} d z = \int_{C^{+}} \sqrt{z} d z + \int_{C^{-}} \sqrt{z} d z = i \int_{ϑ_{b}^{+}}^{ϑ_{f}} e^{i 3 ϑ / 2} d ϑ + i \int_{ϑ_{ι}}^{ϑ_{b}^{-}} \\ e^{i 3 ϑ / 2} d ϑ = \frac{2}{3} [e^{i 3 ϑ_{f} / 2} - e^{i 3 ϑ_{ι} / 2} + e^{i 3 ϑ_{b}^{-} / 2} - e^{i 3 ϑ_{b}^{+} / 2}] = - \frac{4 i}{3}, \end{gathered}

(A68c)

\oint_{C^{'}} \sqrt{z} d z = \frac{2}{3} [e^{i 3 ϑ_{b}^{-} / 2} - e^{i 3 ϑ_{b}^{+} / 2}] = - \frac{4 i}{3},

(A68d)

where we have used that dz = ie^iϑdϑ. It should be underlined that arc C is not well-defined parametrically here, owing to the discontinuity of polar angle ϑ around the negative real axis. In this vein, it contradicts the typical definitions of topological curves, according to which such geometrical entities are continuous mappings from an interval

I

of real numbers into a topological space

X

⁠, i.e.,

C : I \to X

{C is defined here not over a single interval but instead over a union of intervals (−π, −2π/3] ∪ [2π/3, π]}. Moreover, the discontinuous “jump” of ϑ (ϑ = π⁻ → −π⁺) when C intersects the axis ϑ = π leads to issues in the definition of the differential quantity dϑ and subsequently of relation dz = ie^iϑdϑ. All the preceding complications can be simply resolved by considering C as made up of curves C⁺ and C⁻, which are appropriately parameterized over the polar angle intervals [2π/3, π] and (−π, −2π/3], correspondingly. The associated line integrals can then be calculated separately along C⁺ and C⁻, with all differentials being well-defined [a similar rationale was applied in Eqs. and , but there, the separate evaluation along C⁺ and C⁻ was dictated by the branch cuts that either the antiderivative or the integrand quantity exhibited at ϑ = π]. Meanwhile, an analogous approach was not required when computing the contour integrals along the closed path C′, as it is properly defined over (−π, π] and leads to well-behaved differentials [in order to be mathematically rigorous and owing to the presence of the branch cut at ϑ = π, we take the limiting behavior of the aforementioned integrals as their lower and upper bounds tend toward values −π

(θ_{b}^{+} = - π^{+})

and π

(θ_{b}^{-} = π^{-})

⁠, respectively].

We shall conclude the present discussion, by providing the following remarks:

In the current analysis, we have considered the impact of using the principal branch when evaluating line integrals involving multivalued complex functions (the presented approach is general enough and can be applied if other branches are employed). Yet, special care must be taken even when performing more basic computations. For instance, the relation $\sqrt{z_{1} z_{2}} = \sqrt{z_{1}} \sqrt{z_{2}}$ is not always applicable and an example of this is if z₁ = −1 and z₂ = i,
$\begin{gathered} \sqrt{z_{1} z_{2}} = \sqrt{- i} = e^{- i π / 4}, \\ \sqrt{z_{1}} \sqrt{z_{2}} = \sqrt{- 1} \sqrt{i} = e^{i π / 2} e^{i π / 4} = e^{i 3 π / 4} . \end{gathered}$
Similarly, equality ln(z₁z₂) = ln(z₁) + ln(z₂) does not uphold, and to see this, we can set again z₁ = −1 and z₂ = i,
$\begin{gathered} \ln (z_{1} z_{2}) = \ln (- i) = - π / 2, \\ \ln (z_{1}) + \ln (z_{2}) = i (π + π / 2) = i 3 π / 2 . \end{gathered}$
More extensive operations of this type (especially when it comes to square and nth roots) can be found in Sec. II B. In general, to avoid any inconsistencies, all calculations must be done step by step to ensure that each time the principal branches of the different multivalued functions are employed.
Clearly, all integrals evaluated here are well-defined, since they result in finite numerical values, which are also independent of the chosen parameterization. However, in general, for the contour integral ∫_Cf(z)dz to be properly defined, it is necessary and sufficient that C is a piecewise smooth curve and f(z) is piecewise continuous along C. Both of the latter conditions are satisfied in all examined cases.
Certain benefits of employing the $\tilde{(\pm)}$ mode notation, which was introduced in Sec. III, can become apparent. To see this, let us consider the evaluation of the square roots contained in the eigenvalue formulas of Eq. (A3) for a two-mode Hamiltonian arrangement [an analogous rationale can apply for n-mode Hamiltonians, where nth (and lower-order) roots emerge in the corresponding eigenvalue expressions]. When the phase of the complex quantity beneath the square root shifts from the second to the third quadrant (π⁻ → −π⁺) or vice versa (−π⁺ → π⁻), a discontinuity is anticipated in the evolution of λ^±(t) within the complex eigenvalue space given that {z = ρe^iϑ, ρ > 0, with ϑ, ϑ_o ∈ (−π, π]},
$\lim_{ϑ \to ϑ_{o}} \sqrt{z} = \lim_{ϑ \to ϑ_{o}} \sqrt{ρ} e^{i ϑ / 2} = \{\begin{matrix} i \sqrt{ρ}, ϑ_{o} = π^{-}, \\ - i \sqrt{ρ}, ϑ_{o} = - π^{+} . \end{matrix}$
Indeed, in the special case of a $P T$ -symmetric-like system where $λ^{\pm} (t) = \pm i \sqrt{κ^{2} - f^{2} (t)}$ ⁠, such an abrupt transition takes place when λ^±(t) cross the real axis as confirmed from Fig. 2. Meanwhile, when using the $\tilde{(\pm)}$ representation, the trajectories that $λ^{\tilde{\pm}} (t)$ follow remain continuous in any circumstance (see Figs. 2 and 4–6) and this is of course attributed to the definition of the $\tilde{(\pm)}$ notation [i.e., it reorders modes (±) at each time step so as to attain a smooth eigenspectra evolution]. The usefulness then of such a mode formalism is multifold: (i) it helps us to better visualize the evolution of both eigenvalues (e.g., see Fig. 2) and eigenvectors (e.g., see Fig. 3 depicting the eigenvector component ratios), as abrupt “jumps” will not take place in the respective diagrams, (ii) it facilitates the analysis at many points, since it leads to well-defined derivatives and to more manageable integrals, where we do not have to worry about the emergence of branch cuts and where the various theorems/lemmas of complex analysis can become directly applicable, (iii) it yields physically consistent results, as the eigenspectra of an actual physical setting can never exhibit a discontinuous behavior (the eigenstates should vary in a continuous fashion, and the same applies for the associated mode population dynamics, or, more specifically, for the projections of the state vector on the system eigenvectors, assuming of course that the state vector also varies continuously).

APPENDIX B: RIGOROUS ANALYTICAL TREATMENT OF THE MODE POPULATION DYNAMICS

1. Deriving the mode population equations for generalized Hamiltonian settings

At this point, we shall retrieve the mode population dynamics for two-level Hamiltonian arrangements. We will be using the scaled time coordinate $\bar{t}$ (=t/T), so as (i) the dependencies of the different variables and parameters from $\bar{t}$ and T can be clearly depicted and (ii) to facilitate the slowly varying study in Appendix B 3, where T serves as a tuning knob to achieve adiabaticity. Moreover, the $(\tilde{\pm})$ mode formalism shall be employed in order to avoid potential discontinuities that may affect the evaluation of complex integrals and, additionally, lead to ill-defined derivatives (see the relevant discussion in Appendix A 5).

To initiate our analysis, let us decompose the state vector

|u (t)⟩

or

|u^{T} (\bar{t})⟩

(for notation details, see Appendix A 2 and particularly the second bullet point toward the end of Appendix A 2) in the biorthonormal basis

{|V_{R}^{\tilde{+}} (\bar{t})⟩, |V_{R}^{\tilde{-}} (\bar{t})⟩}

as follows:

\begin{align} |u^{T} (\bar{t})⟩ & = c_{T}^{\tilde{+}} (\bar{t}) e^{T \int_{0}^{\bar{t}} λ^{\tilde{+}} ({\bar{t}}^{'}) d {\bar{t}}^{'}} |V_{R}^{\tilde{+}} (\bar{t})⟩ + c_{T}^{\tilde{-}} (\bar{t}) e^{T \int_{0}^{\bar{t}} λ^{\tilde{-}} ({\bar{t}}^{'}) d {\bar{t}}^{'}} |V_{R}^{\tilde{-}} (\bar{t})⟩, \end{align}

(B1)

where integration from 0 to

\bar{t}

has been assumed. [The preceding expression is equivalent to Eq. in Sec. III, keeping in mind equality

\int_{0}^{t} λ^{\tilde{\pm}} (t^{'}) d t^{'} = T \int_{0}^{\bar{t}} λ^{\tilde{\pm}} ({\bar{t}}^{'}) d {\bar{t}}^{'}

⁠.] Substituting Eq. in Eq.