Quasiclassical analysis of Bloch oscillations in non-Hermitian tight-binding lattices

Non-Hermitian quantum mechanics has in recent years generated substantial amounts of research interest. In particular the realisation of PT-symmetric quantum dynamics in the context of optics [1–3] has opened up a whole new field of investigations. In this context, periodic potentials play an important role and have been studied extensively in the literature [4–10]. One of the most striking features of unitary quantum dynamics in periodic potentials is the phenomenon of Bloch oscillations [11–13]. If a static force is applied to a periodic lattice, instead of being transported in the direction of the force, quantum particles perform the famous Bloch oscillations, which have been observed in a variety of experimental systems reaching from semiconductor superlattices to optical waveguide structures [14–18]. Recently the effect of static forces and the modification of Bloch oscillations in PT-symmetric and more general non-Hermitian lattices has been investigated both theoretically and experimentally [6, 9, 19, 20].

Interestingly, in the unitary case, despite their counterintiuitive nature, Bloch oscillations can be understood on the basis of a simple quasiclassical argument building on Hamilton's equations of motion. Many features of the exact quantum dynamics can be qualitatively recovered from such a quasiclassical description. The dynamics of the classical counterpart of non-Hermitian quantum systems, however, is more subtle. It has been recently argued in [9] that a quasiclassical description of Bloch oscillations in the non-Hermitian case is of little use. Here, however, we show that the classical dynamics recently developed in [21–24] as a counterpart of non-Hermitian quantum dynamics, is capable of describing the main features of Bloch oscillations as well in the non-Hermitian case, as long as only a single Bloch band is involved in the dynamics (a constraint that also holds in the Hermitian case).

The paper is organised as follows: we first introduce the model system studied here, a non-Hermitian single-band tight-binding Hamiltonian, and summarise some of the important properties of the quantum description. We then review the quasiclassical description of Bloch oscillations in the Hermitian case, and the equations of motion arising as a semiclassical limit of quantum dynamics generated by non-Hermitian Hamiltonians. Finally we study two example systems in detail, a Hatano–Nelson lattice with real coupling constants and a system with purely imaginary couplings. It is shown how basic features of the dynamics can be understood from the quasiclassical description both for initially localised states in space or in momentum.

Here we study a non-Hermitian extension of a single-band tight-binding system with Hamiltonian

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{n=-\infty }^{+\infty }{g}_{1}| n\rangle \langle n+1| +{g}_{2}| n+1\rangle \langle n| +2{Fn}| n\rangle \langle n| ,\end{eqnarray} \tag{ 1 }$$

with an orthogonal basis $| n\rangle$ (for example the Wannier states of a periodic potential). The parameters g₁ and g₂ describe the hopping between neighbouring sites and $2F$ is the static force. Here we restrict the discussion to real force parameters $F\in {\mathbb{R}}$. For ${g}_{1}^{* }\;=\;{g}_{2}$ the system is Hermitian, and performs Bloch oscillations (see, e.g., [13] and references therein for more information). In what follows we shall use units with ${\hslash }=1$, and confine ourselves to time-independent parameters. While in the following we adopt the language of quantum mechanics, the Hamiltonian (1) can be realised in an optical setup using waveguides or fibre loops (see [9, 10] and references therein).

Let us start by reviewing some elementary properties of the system described by the Hamiltonian (1) based on the well-known results for the Hermitian case (see, e.g., [13] and references given there). It is convenient to define the operators

$$\begin{eqnarray}&&\hat{K}=\displaystyle \sum _{n}| n\rangle \langle n+1| ,\quad {\hat{K}}^{\dagger }=\displaystyle \sum _{n}| n+1\rangle \langle n| ,\quad {\rm{and}}\quad \hat{N}=\displaystyle \sum _{n}n| n\rangle \langle n| \end{eqnarray} \tag{ 2 }$$

with ${\hat{N}}^{\dagger }=\hat{N}$. The operators $\hat{K}$ and ${\hat{K}}^{\dagger }$ are unitary and act on the $| n\rangle$ states as ladder or shift operators, i.e.,

$$\begin{eqnarray}&&\hat{K}| n\rangle =| n-1\rangle ,\quad {\hat{K}}^{\dagger }| n\rangle =| n+1\rangle .\end{eqnarray} \tag{ 3 }$$

The operators (2) satisfy the commutation relations

$$\begin{eqnarray}&&[\hat{K},\hat{N}]=\hat{K},\quad [{\hat{K}}^{\dagger },\hat{N}]=-{\hat{K}}^{\dagger },\quad [\hat{K},{\hat{K}}^{\dagger }]=0,\end{eqnarray} \tag{ 4 }$$

and form a Lie algebra denoted as the shift operator algebra [25, 26]. The Hamiltonian (1) written in terms of the operators (2) reads

$$\begin{eqnarray}&&\hat{H}={g}_{1}\hat{K}+{g}_{2}{\hat{K}}^{\dagger }+2F\hat{N}.\end{eqnarray} \tag{ 5 }$$

It is often convenient to introduce the quasi-momentum operator $\hat{\kappa }$ defined via

$$\begin{eqnarray}&&\hat{K}={{\rm{e}}}^{{\rm{i}}\hat{\kappa }},\end{eqnarray} \tag{ 6 }$$

as the conjugate observable of the discrete position operator $\hat{N}$. They fulfil the canonical commutation relation

$$\begin{eqnarray}&&[\hat{N},\hat{\kappa }]={\rm{i}}.\end{eqnarray} \tag{ 7 }$$

In these operators the Hamiltonian can be written as

$$\begin{eqnarray}&&\hat{H}=E(\hat{\kappa })+2F\hat{N},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&E(\hat{\kappa })={g}_{1}{{\rm{e}}}^{{\rm{i}}\hat{\kappa }}+{g}_{2}{{\rm{e}}}^{-{\rm{i}}\hat{\kappa }}\end{eqnarray} \tag{ 9 }$$

is the dispersion relation for the field free case F = 0. In what follows it is convenient to decompose the Hamiltonian into Hermitian and anti-Hermitian parts as $\hat{H}={\hat{H}}_{R}-{\rm{i}}{\hat{H}}_{I}$, with

$$\begin{eqnarray}&&{\hat{H}}_{R}=\mathrm{Re}\;{g}_{+}\mathrm{cos}(\hat{\kappa })-\mathrm{Im}\;{g}_{-}\mathrm{sin}(\hat{\kappa })+2F\hat{N},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\hat{H}}_{I}=-\mathrm{Im}\;{g}_{+}\mathrm{cos}(\hat{\kappa })-\mathrm{Re}\;{g}_{-}\mathrm{sin}(\hat{\kappa }),\end{eqnarray} \tag{ 11 }$$

where we have defined ${g}_{\pm }\;=\;{g}_{1}\pm {g}_{2}$. It is straightforward to verify that for real values of g_1,2 the Hamiltonian is PT-symmetric, i.e., invariant under the PT-transformation $\kappa \to -\kappa$, ${\rm{i}}\to -{\rm{i}}$, and $N\to N$ ³ . The Hamiltonian is in fact PT-symmetric under the more general transformation $\kappa \to \theta -\kappa$, ${\rm{i}}\to -{\rm{i}}$, $N\to N$, whenever the two complex parameters g_1,2 have the same phase θ.

The Bloch waves

$$\begin{eqnarray}&&| \kappa \rangle =\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \sum _{n}{{\rm{e}}}^{{\rm{i}}n\kappa }| n\rangle \end{eqnarray} \tag{ 12 }$$

are eigenstates of $\hat{K}$ and ${\hat{K}}^{\dagger }$ with

$$\begin{eqnarray}&&\hat{K}| \kappa \rangle ={{\rm{e}}}^{{\rm{i}}\kappa }| \kappa \rangle ,\quad {\hat{K}}^{\dagger }| \kappa \rangle ={{\rm{e}}}^{-{\rm{i}}\kappa }| \kappa \rangle ,\end{eqnarray} \tag{ 13 }$$

where κ is the quasimomentum. They are orthogonal and normalised to $2\pi$-periodic delta-comb functions as

$$\begin{eqnarray}&&\langle \kappa | \kappa ^{\prime} \rangle =\displaystyle \frac{1}{2\pi }\displaystyle \sum _{n}{{\rm{e}}}^{{\rm{i}}n(\kappa ^{\prime} -\kappa )}={\delta }_{2\pi }(\kappa ^{\prime} -\kappa ).\end{eqnarray} \tag{ 14 }$$

The matrix elements of the position operator in the Bloch basis are given by

$$\begin{eqnarray}&&\langle \kappa | \hat{N}| \kappa ^{\prime} \rangle =-{\rm{i}}{\delta }_{2\pi }^{\prime }(\kappa ^{\prime} -\kappa )={\rm{i}}{\delta }_{2\pi }(\kappa ^{\prime} -\kappa )\displaystyle \frac{{\rm{d}}\ }{{\rm{d}}\kappa }.\end{eqnarray} \tag{ 15 }$$

Therefore, the Hamiltonian is diagonal in the quasimomentum representation, $\langle \kappa | \hat{H}| \kappa ^{\prime} \rangle =\hat{H}(\kappa )\delta (\kappa ^{\prime} -\kappa )$ with

$$\begin{eqnarray}&&\hat{H}(\kappa )={g}_{1}{{\rm{e}}}^{{\rm{i}}\kappa }+{g}_{2}{{\rm{e}}}^{-{\rm{i}}\kappa }+{\rm{i}}\;2F\displaystyle \frac{{\rm{d}}}{{\rm{d}}\kappa }.\end{eqnarray} \tag{ 16 }$$

For time-independent parameters, the eigenstates of this Hamiltonian for periodic boundary conditions ${\rm{\Psi }}(0)={\rm{\Psi }}(2\pi )$ are

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{m}(\kappa )\propto {{\rm{e}}}^{-{\rm{i}}m\kappa +\tfrac{{g}_{1}}{2F}{{\rm{e}}}^{{\rm{i}}\kappa }-\displaystyle \frac{{g}_{2}}{2F}{{\rm{e}}}^{-{\rm{i}}\kappa }}\end{eqnarray} \tag{ 17 }$$

with eigenvalues

$$\begin{eqnarray}&&{E}_{m}=2{Fm},\quad m=0,\pm 1,\pm 2,\ldots .\end{eqnarray} \tag{ 18 }$$

It should be noted that these eigenvalues are the same as for the Hermitian case, they are real valued, equidistant and independent of the coefficients g_j as has already been observed in [9]. Due to the equidistance of the eigenvalues the quantum dynamics is necessarily periodic with Bloch period ${T}_{{\rm{B}}}=\pi /F$. Nevertheless the dynamics is influenced by the non-Hermiticity in a nontrivial manner, due to the modification of eigenstates (17), and in particular due to the fact that they are not orthogonal in general.

In terms of the coefficients c_n of the wave function in the basis $| n\rangle$ the time-dependent Schrödinger equation reads

$$\begin{eqnarray}&&{\rm{i}}\;\displaystyle \frac{{\rm{d}}{c}_{n}}{{\rm{d}}t}={g}_{1}{c}_{n+1}+{g}_{2}{c}_{n-1}+2{{Fnc}}_{n},\end{eqnarray} \tag{ 19 }$$

which can of course be solved numerically. In the Hermitian case (i.e., ${g}_{2}={g}_{1}^{* }$), the time-evolution can be obtained analytically using algebraic techniques [13]. We note that these techniques can be extended to the non-Hermitian case, and analytic results for the time-evolution can be obtained. Here, however, we focus on the quasiclassical description of the dynamics.

In the Hermitian case the Heisenberg equations of motion for the expectation values of the quasimomentum and discrete position operators are given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle \hat{\kappa }\rangle =-2F\quad {\rm{and}}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle \hat{N}\rangle =\left\langle \displaystyle \frac{\partial E(\hat{\kappa })}{\partial \hat{\kappa }}\right\rangle ,\end{eqnarray} \tag{ 21 }$$

where $E(\hat{\kappa })$ is the field free dispersion relation (9). The equation of motion for $\langle \hat{\kappa }\rangle$ is trivially integrated to yield a linearly changing quasimomentum, $\langle \hat{\kappa }\rangle (t)=-2{Ft}+\langle \hat{\kappa }\rangle (0)$, a result known as the acceleration theorem [27]. In the spirit of Ehrenfest's theorem one can replace the (real valued) expectation value $\left\langle \tfrac{\partial E(\hat{\kappa })}{\partial \hat{\kappa }}\right\rangle$ with $\tfrac{E(\langle \hat{\kappa }\rangle )}{\partial \langle \hat{\kappa }\rangle }$ to obtain the quasiclassical equations of motion given by Hamilton's equations of motion

$$\begin{eqnarray}&&\dot{p}=-\displaystyle \frac{\partial H}{\partial q}=-2F\quad {\rm{and}}\quad \dot{q}=\displaystyle \frac{\partial H}{\partial p}=\displaystyle \frac{\partial E(p)}{\partial p},\end{eqnarray} \tag{ 22 }$$

where we have introduced the classical momentum and position variables p and q as the expectation values of $\hat{\kappa }$ and $\hat{N}$ respectively, and with the classical Hamiltonian function

$$\begin{eqnarray}&&H={g}_{1}{{\rm{e}}}^{{\rm{i}}p}+{g}_{2}{{\rm{e}}}^{-{\rm{i}}p}+2{Fq},\end{eqnarray} \tag{ 23 }$$

with ${g}_{2}={g}_{1}^{* }$. The equation for the position variable is integrated to yield

$$\begin{eqnarray}&&q(t)={q}_{0}+\displaystyle \frac{E({p}_{0})-E(p(t))}{2F},\end{eqnarray} \tag{ 24 }$$

with $p(t)={p}_{0}-2{Ft}$, where ${p}_{0}=p(0)$ and equivalently for q. That is, the position variable maps out the periodic dispersion relation over time—which is the quasiclassical explanation of the phenomenon of Bloch oscillations [9, 13, 28].

In the non-Hermitian case, the equation of motion for expectation values $\langle \hat{A}\rangle =\tfrac{\langle \psi | \hat{A}| \psi \rangle }{\langle \psi | \psi \rangle }$ is modified in a non-trivial way to [21, 29, 30]

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle \hat{A}\rangle =\langle [\hat{A},{\hat{H}}_{R}]\rangle -{\rm{i}}(\langle {[\hat{A},{\hat{H}}_{I}]}_{+}\rangle -2\langle \hat{A}\rangle \langle {\hat{H}}_{I}\rangle ).\end{eqnarray} \tag{ 25 }$$

That is, the dynamical equations for observables depend on the covariances with the anti-Hermitian part of the Hamiltonian. In [9] it has been argued that the resulting dynamics for the variables p and q are of little use for the understanding of the dynamics in the non-Hermitian case. In the following we shall show, however, that while it may indeed be difficult to proceed directly from equation (25) for position and momentum, the semiclassical limit of the quantum dynamics generated by non-Hermitian Hamiltonians, as derived in [22] is capable of accurately describing basic features of the quantum dynamics.

Using a Gaussian wavepacket approximation in the spirit of Heller [31], it has been shown that the classical dynamics associated to a non-Hermitian Hamiltonian $H={H}_{R}-{\rm{i}}{H}_{I}$, depending on the real valued canonical coordinates $p,q$ are given by [22, 23]

$$\begin{eqnarray}&&\dot{p}=-\displaystyle \frac{\partial {H}_{R}}{\partial q}-{{\rm{\Sigma }}}_{{pp}}\displaystyle \frac{\partial {H}_{I}}{\partial p}-{{\rm{\Sigma }}}_{{pq}}\displaystyle \frac{\partial {H}_{I}}{\partial q},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\dot{q}=\displaystyle \frac{\partial {H}_{R}}{\partial p}-{{\rm{\Sigma }}}_{{pq}}\displaystyle \frac{\partial {H}_{I}}{\partial p}-{{\rm{\Sigma }}}_{{qq}}\displaystyle \frac{\partial {H}_{I}}{\partial q},\end{eqnarray} \tag{ 27 }$$

where Σ is a real symmetric matrix, proportional to the covariance matrix, that is, it encodes the (co)variances of position and momentum according to

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{pp}}=\displaystyle \frac{2}{{\hslash }}{({\rm{\Delta }}p)}^{2},\quad {{\rm{\Sigma }}}_{{qq}}=\displaystyle \frac{2}{{\hslash }}{({\rm{\Delta }}q)}^{2},\quad {\rm{and}}\quad {{\rm{\Sigma }}}_{{pq}}={{\rm{\Sigma }}}_{{qp}}=\displaystyle \frac{2}{{\hslash }}{{\rm{\Delta }}}_{{pq}},\end{eqnarray} \tag{ 28 }$$

where the determinant of Σ is one. Note that we have included the ℏ here for clarity, while we shall continue to use rescaled units with ${\hslash }=1$ in the following. The covariances are also time-dependent, following the classical dynamical equations

$$\begin{eqnarray}&&\dot{{\rm{\Sigma }}}={\rm{\Omega }}{H}_{R}^{\prime\prime }{\rm{\Sigma }}-{\rm{\Sigma }}{H}_{R}^{\prime\prime }{\rm{\Omega }}-{\rm{\Omega }}{H}_{I}^{\prime\prime }{\rm{\Omega }}-{\rm{\Sigma }}{H}_{I}^{\prime\prime }{\rm{\Sigma }},\end{eqnarray} \tag{ 29 }$$

where Ω is the standard symplectic matrix

$$\begin{eqnarray}{\rm{\Omega }}=\left(\begin{array}{ll}0 & -1\\ 1 & 0\end{array}\right),\end{eqnarray} \tag{ 30 }$$

and ${H}_{R}^{\prime\prime }$ and ${H}_{I}^{\prime\prime }$ denote the matrices of second phase-space derivatives of H_R and H_I respectively. The most striking difference to standard Hamiltonian dynamics is the dissipative term in the equations of motion for the canonical variables that couples to the dynamics of the width parameters. In addition we have a classical approximation for the time evolution of the squared norm (or the total power, in the context of optics) $P=\langle \psi | \psi \rangle$ of the wave packet given by

$$\begin{eqnarray}&&\dot{P}=-\left(2{H}_{I}-\displaystyle \frac{1}{2}\mathrm{Tr}({\rm{\Omega }}{H}_{I}^{\prime\prime }{\rm{\Omega }}{{\rm{\Sigma }}}^{-1})\right)P.\end{eqnarray} \tag{ 31 }$$

Note that we constrain the discussion to one-dimensional systems here, and use the scaled covariance matrix Σ instead of its inverse G that is used in [22, 23].

In the present case for general parameter values we have

$$\begin{eqnarray}&&{H}_{R}=\mathrm{Re}\;{g}_{+}\mathrm{cos}\;p-\mathrm{Im}\;{g}_{-}\mathrm{sin}\;p+2{Fq},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{H}_{I}=-\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p-\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p.\end{eqnarray} \tag{ 33 }$$

Thus, most of the elements of the Hessian matrices ${H}_{I,R}^{\prime\prime }$ are zero, and we find the relatively compact classical equations of motion

$$\begin{eqnarray}&&\dot{p}=-2F-(\mathrm{Im}\;{g}_{+}\mathrm{sin}\;p-\mathrm{Re}\;{g}_{-}\mathrm{cos}\;p){{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\dot{q}=-\mathrm{Re}\;{g}_{+}\mathrm{sin}\;p-\mathrm{Im}\;{g}_{-}\mathrm{cos}\;p-(\mathrm{Im}\;{g}_{+}\mathrm{sin}\;p-\mathrm{Re}\;{g}_{-}\mathrm{cos}\;p){{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pp}}=-(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p){{\rm{\Sigma }}}_{{pp}}^{2},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pq}}=(-\mathrm{Re}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Im}\;{g}_{-}\mathrm{sin}\;p-(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p){{\rm{\Sigma }}}_{{pq}}){{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{qq}}=-2(\mathrm{Re}\;{g}_{+}\mathrm{cos}\;p-\mathrm{Im}\;{g}_{-}\mathrm{sin}\;p){{\rm{\Sigma }}}_{{pq}}+(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p)(1-{{\rm{\Sigma }}}_{{pq}}^{2}),\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&\dot{P}=(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p)\left(2-\tfrac{1}{2}{{\rm{\Sigma }}}_{{pp}}\right)P,\end{eqnarray} \tag{ 39 }$$

where again we have used the notation ${g}_{\pm }\;=\;{g}_{1}\pm {g}_{2}$.

The classical dynamics are expected to be a good approximation of the quantum dynamics as long as the width in momentum space (described by ${{\rm{\Sigma }}}_{{pp}}$) stays small, since the Hamiltonian is anharmonic in momentum. It is interesting to note that indeed if ${{\rm{\Sigma }}}_{{pp}}$ is negligible initially, according to equation (36) it will stay negligible during the time evolution. In this case we can perform a further approximation, by setting ${{\rm{\Sigma }}}_{{pp}}$ to zero in the equations of motion, from which it follows that ${{\rm{\Sigma }}}_{{pq}}$ stays constant. Thus, we find the approximate dynamics for narrow momentum wave packets given by

$$\begin{eqnarray}&&\dot{p}=-2F,\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\dot{q}=-\mathrm{Re}\;{g}_{+}\mathrm{sin}\;p-\mathrm{Im}\;{g}_{-}\mathrm{cos}\;p-(\mathrm{Im}\;{g}_{+}\mathrm{sin}\;p-\mathrm{Re}\;{g}_{-}\mathrm{cos}\;p){{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{qq}}=-2(\mathrm{Re}\;{g}_{+}\mathrm{cos}\;p-\mathrm{Im}\;{g}_{-}\mathrm{sin}\;p){{\rm{\Sigma }}}_{{pq}}+(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p)(1-{{\rm{\Sigma }}}_{{pq}}^{2}),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\dot{P}=\ 2(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p)P.\end{eqnarray} \tag{ 43 }$$

In this approximation the acceleration theorem still holds and the momentum changes linearly with time according to

$$\begin{eqnarray}&&p(t)={p}_{0}-2{Ft}.\end{eqnarray} \tag{ 44 }$$

The dynamical equation for the position can be rewritten as

$$\begin{eqnarray}&&\dot{q}=\displaystyle \frac{\partial \mathrm{Re}\;(E(p))}{\partial p}+\displaystyle \frac{\partial \mathrm{Im}\;(E(p))}{\partial p}{{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 45 }$$

where E(p) is the field-free dispersion relation given by

$$\begin{eqnarray}E(p) & = & {g}_{1}{{\rm{e}}}^{{\rm{i}}p}+{g}_{2}{{\rm{e}}}^{-{\rm{i}}p}\\ & = & \mathrm{Re}\;{g}_{+}\mathrm{cos}\;p-\mathrm{Im}\;{g}_{-}\mathrm{sin}\;p+{\rm{i}}(\mathrm{Im}\;{g}_{+}\mathrm{cos}\;p+\mathrm{Re}\;{g}_{-}\mathrm{sin}\;p),\end{eqnarray} \tag{ 46 }$$

and where the covariance ${{\rm{\Sigma }}}_{{pq}}$ stays constant in time. That is, for vanishing covariance the motion in space follows the real part of the field-free dispersion relation, similar to the Hermitian case, as has been observed in [9]. For non-vanishing covariance of p and q, however, the imaginary part of the dispersion relation also contributes to the dynamics. In general the time dependence of the position for narrow momentum wave packets is given by

$$\begin{eqnarray}&&q(t)={q}_{0}-\tfrac{\mathrm{Re}\;{g}_{+}+\mathrm{Im}\;{g}_{+}{{\rm{\Sigma }}}_{{pq}}}{2F}(\mathrm{cos}\;p(t)-\mathrm{cos}\;{p}_{0})+\tfrac{\mathrm{Im}\;{g}_{-}-\mathrm{Re}\;{g}_{-}{{\rm{\Sigma }}}_{{pq}}}{2F}(\mathrm{sin}\;p(t)-\mathrm{sin}\;{p}_{0}),\end{eqnarray} \tag{ 47 }$$

and the squared norm evolves as

$$\begin{eqnarray}&&P(t)=\mathrm{exp}\left(-\tfrac{\mathrm{Im}\;{g}_{+}}{F}(\mathrm{sin}\;p(t)-\mathrm{sin}\;{p}_{0})-\tfrac{\mathrm{Re}\;{g}_{-}}{F}(\mathrm{cos}\;p(t)-\mathrm{cos}\;{p}_{0})\right){P}_{0}.\end{eqnarray} \tag{ 48 }$$

For initial states with small but non zero width in momentum, the classical dynamics is modified from this simple picture. Note that the classical equations in general do not lead to periodic motions in space, in contrast to the full quantum dynamics.

It may at first seem impossible to use the quasiclassical description for wave packets that are narrow in position space, since a small value of ${{\rm{\Sigma }}}_{{qq}}$ implies a large value of ${{\rm{\Sigma }}}_{{pp}}$ due to the uncertainty relation. However, we can decompose a narrow wave packet in position space into wave packets that are narrow in momentum space and broad in position. In particular, we can interpret the Fourier transform of a narrow wave packet in position space as a superposition of plane waves, which we can then approximate by a classical ensemble of plane waves with the appropriate momentum distribution. Each of these plane waves is an extreme case of a Gaussian wave packet with ${{\rm{\Sigma }}}_{{pp}}=0$ centred at q = 0 and arbitrary p. Following this argument each of the individual ensemble trajectories evolves in time according to equations (40)–(42). They must be weighted with their respective squared norm which evolves according to equation (43). That is, we replace the weighted average of a variable A for a state localised initially in position q, denoted by $\langle A{\rangle }^{(q)}$, by the integral over the weighted average of states localised initially in momentum, denoted by $\langle A{\rangle }^{(p)}$, that is

$$\begin{eqnarray}&&{P}^{(q)}\langle A{\rangle }^{(q)}=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{\rm{d}}p\;{P}^{(p)}\langle A{\rangle }^{(p)}.\end{eqnarray} \tag{ 49 }$$

This argument provides a remarkably good approximation of the dynamics and in fact reproduces the exact quantum results for the mean value of the position and the squared norm of the wave packet in the case of initial wave packet localised in a single site, as will be demonstrated for the two examples discussed in the following. In an appendix we provide an analytical quantum mechanical argument why the ensemble method is expected to yield exact results for the expectation values of operators that are within the algebra spanned by $\hat{K},{\hat{K}}^{\dagger },\hat{N}$, and the identity operator.

To demonstrate the strength of the quasiclassical description we consider two PT-symmetric examples with non-Hermitian coupling constants leading to complex dispersion relations that have been discussed in [9].

4.1. Hatano–Nelson lattice

As a first example we consider a model inspired by the Hatano–Nelson Hamiltonian [32], which was introduced in the context of magnetic flux lines in superconductors. In this model the coupling constants are real, but asymmetric, ${g}_{\mathrm{1,2}}=g{{\rm{e}}}^{\pm \mu }$, with $g,\mu \in {\mathbb{R}}$ describing a biased coupling. That is, the Hamiltonian is given by

$$\begin{eqnarray}&&H=2g\mathrm{cosh}\mu \mathrm{cos}\;p+2{\rm{i}}g\mathrm{sinh}\mu \mathrm{sin}\;p+2{Fq}.\end{eqnarray} \tag{ 50 }$$

An experimental implementation of this model using optical resonator structures has been proposed in [10]. The quantum model can be mapped to an isospectral Hermitian Hamiltonian by the substitution ${c}_{n}={{\rm{e}}}^{-\mu n}{\tilde{c}}_{n}$ with [9]

$$\begin{eqnarray}&&{\rm{i}}\;\displaystyle \frac{{\rm{d}}{\tilde{c}}_{n}}{{\rm{d}}t}=g\;{\tilde{c}}_{n+1}+g\;{\tilde{c}}_{n-1}+2{Fn}{\tilde{c}}_{n}.\end{eqnarray} \tag{ 51 }$$

Therefore, the solution of the quantum dynamics can be obtained from that of the Hermitian system (51), the time evolution matrix of which is given by

$$\begin{eqnarray}&&{\tilde{U}}_{{nn}^{\prime} }(t)={J}_{n-n^{\prime} }\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{{\rm{i}}(n-n^{\prime} )(\tfrac{\pi }{2}-{Ft})-{\rm{i}}n^{\prime} 2{Ft}},\end{eqnarray} \tag{ 52 }$$

where J_m(z) is the Bessel function. Hence the time evolution matrix for the non-Hermitian case is simply given by [9]

$$\begin{eqnarray}{U}_{{nn}^{\prime} }(t) & = & {{\rm{e}}}^{-\mu (n-n^{\prime} )}{\tilde{U}}_{{nn}^{\prime} }(t)\\ & = & {J}_{n-n^{\prime} }\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{{\rm{i}}(n-n^{\prime} )(\tfrac{\pi }{2}-{Ft}+{\rm{i}}\mu )-{\rm{i}}n^{\prime} 2{Ft}}\end{eqnarray} \tag{ 53 }$$

and the c_n(t) coefficients are

$$\begin{eqnarray}&&{c}_{n}(t)=\displaystyle \sum _{n^{\prime} }{U}_{{nn}^{\prime} }(t){c}_{n}(0).\end{eqnarray} \tag{ 54 }$$

The classical equations of motion are given by

$$\begin{eqnarray}&&\dot{p}=-2F+2g\mathrm{sinh}\mu \mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&\dot{q}=-2g\mathrm{cosh}\mu \mathrm{sin}\;p+2g\mathrm{sinh}\mu \mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pp}}=-2g\mathrm{sinh}\mu \mathrm{sin}\;p\;{{\rm{\Sigma }}}_{{pp}}^{2},\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pq}}=-2g(\mathrm{cosh}\mu \mathrm{cos}\;p+\mathrm{sinh}\mu \mathrm{sin}\;p\;{{\rm{\Sigma }}}_{{pq}}){{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{qq}}=2g\mathrm{sinh}\mu \mathrm{sin}\;p(1-{{\rm{\Sigma }}}_{{pq}}^{2})-4g\mathrm{cosh}\mu \mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&\dot{P}=-g\mathrm{sinh}\mu \mathrm{sin}\;p({{\rm{\Sigma }}}_{{pp}}-4)P.\end{eqnarray} \tag{ 60 }$$

In the limiting case of infinitely narrow momentum wave packets with ${{\rm{\Sigma }}}_{{pp}}=0$ the momentum moves according to the acceleration theorem (44), and we find the simplified dynamical equations

$$\begin{eqnarray}&&\dot{q}=-2g\mathrm{cosh}\mu \mathrm{sin}\;p+2g\mathrm{sinh}\mu \mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\dot{P}=4g\mathrm{sinh}\mu \mathrm{sin}\;p\;P,\end{eqnarray} \tag{ 62 }$$

for the position and the squared norm. That is, q performs oscillations

$$\begin{eqnarray}&&q(t)={q}_{0}-\tfrac{2g\mathrm{cosh}\mu }{F}\mathrm{sin}({Ft})\mathrm{sin}({p}_{0}-{Ft})-\tfrac{2g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\mathrm{cos}({p}_{0}-{Ft}){{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 63 }$$

which for vanishing initial covariance ${{\rm{\Sigma }}}_{{pq}}$ between position and momentum follow the real part of the dispersion relation, and are only modified by a scaling factor depending on the non-Hermitian parameter μ in comparison to the Hermitian case $\mu =0$. For non-vanishing covariance there is a qualitative change in the dynamics compared to the Hermitian case. The squared norm is in both cases strongly influenced by the non-Hermiticity and varies according to

$$\begin{eqnarray}&&P(t)=\mathrm{exp}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\mathrm{sin}({p}_{0}-{Ft})\right){P}_{0}.\end{eqnarray} \tag{ 64 }$$

This simple picture gives a good understanding of the quantum dynamics of broad wave packets, as demonstrated in figure 1. In the left panel we show the (renormalised) quantum evolution of an initial broad Gaussian wave packet

$$\begin{eqnarray}&&{c}_{n}(0)\propto {{\rm{e}}}^{-\beta {n}^{2}},\end{eqnarray} \tag{ 65 }$$

where the parameter β encodes the covariances of the wave packet according to

$$\begin{eqnarray}&&\beta =\displaystyle \frac{1}{2{{\rm{\Sigma }}}_{{qq}}}(1-{\rm{i}}{{\rm{\Sigma }}}_{{pq}}),\end{eqnarray} \tag{ 66 }$$

and which we choose as $\beta =0.02$ in the top and $\beta =0.004(1-2{\rm{i}})$ in the bottom panel, corresponding to a broad wave packet with ${q}_{0}=0$, ${p}_{0}=0$, and vanishing and non-vanishing covariance between position and momentum, respectively, with the same small value of ${{\rm{\Sigma }}}_{{pp}}$. The wave packet with vanishing initial covariance performs typical Bloch oscillations similar to the Hermitian case. The packet with non-vanishing initial covariance depicted in the lower panel, also oscillates, but it can be clearly seen that the oscillation is also influenced by the imaginary part of the dispersion relation, leading it to also oscillate to the left. In addition the covariance induces pronounced modulations of the width ${{\rm{\Sigma }}}_{{qq}}$. Note that we have chosen a larger value of μ for the case with non-vanishing covariance, to make the effect more obvious. The middle panel of the figure shows the corresponding classical approximations, in which the wave packet is assumed to stay a Gaussian for all times, described by the distribution

$$\begin{eqnarray}&&{c}_{n}(t)\propto \sqrt{P}{{\rm{e}}}^{-\beta {n}^{2}-{\rm{i}}p(n-q)},\end{eqnarray} \tag{ 67 }$$

where β is related to the covariances via (66), and where the parameters $q,p,{\rm{\Sigma }}$, and P follow the classical dynamics (55)–(60). The initial value of the momentum width described by ${{\rm{\Sigma }}}_{{pp}}$ follows from the condition $\mathrm{det}{\rm{\Sigma }}=1$. The right panel shows the comparison of the motion of the wave packet centre in the full quantum dynamics, the classical approximation, and the simplified classical approximation for infinitely narrow momentum distribution. We observe a good agreement between the quantum and classical descriptions, and both can be well understood in terms of the simplified dynamics for vanishing ${{\rm{\Sigma }}}_{{pp}}$.

Zoom In Zoom Out Reset image size

In figure 2 we depict the squared norm of the propagated beam for two examples, one corresponding to the example on the left in figure 1, and the other to the same initial conditions where the sign of μ is reversed. We note a striking difference in the behaviour for the different signs of μ, belonging to a biased transport to the left or right, respectively. It might be surprising that the evolution of the renormalised beam in the lattice differs only very little for the two cases. This behaviour is in fact well described by the classical approximation, which is shown in comparison to the full quantum evolution in the figure, and agrees perfectly on the depicted scale. We also plot the results obtained from the simplified classical approximation for ${{\rm{\Sigma }}}_{{pp}}$ given by equation (64), in which the opposite behaviour depending on the sign of μ is immediately obvious, and which agrees well with the two more accurate descriptions. In the simplified classical dynamics it is also obvious why the path of the beam in the lattice is not influenced by the sign of μ, since for vanishing initial covariance, it simply maps out the real part of the dispersion relation, which depends only on the absolute value of μ.

Zoom In Zoom Out Reset image size

While much of the dynamics of broad wave packets is well explained using the simplified picture of zero width in momentum, there are visible deviations from the acceleration theorem in the dynamics of p even for the examples shown in figure 1. As has already been observed in [9] the expectation value of the momentum shows small periodic fluctuations around the value ${p}_{0}-2{Ft}$. In figure 3 on the left we show these deviations for the example in figure 1, using both the quantum and the classical descriptions, which agree in this case. We can obtain an analytical approximation for the periodic deviations from the acceleration theorem for small but non-zero ${{\rm{\Sigma }}}_{{pp}}$ from equations (55) and (58) in the following way. Formally integrating equation (55) we find

$$\begin{eqnarray}&&p(t)={p}_{0}-2{Ft}+2g\mathrm{sinh}\mu {\int }_{0}^{t}{{\rm{\Sigma }}}_{{pp}}\mathrm{cos}(p(t)){\rm{d}}t.\end{eqnarray} \tag{ 68 }$$

Integration by parts and using the fact that ${\dot{{\rm{\Sigma }}}}_{{pp}}\sim { \mathcal O }({{\rm{\Sigma }}}_{{pp}}^{2})$ yields

$$\begin{eqnarray}&&p={p}_{0}-2{Ft}+2g\mathrm{sinh}\mu {\left[{{\rm{\Sigma }}}_{{pp}}(\tau )\int \mathrm{cos}(p(\tau )){\rm{d}}\tau \right]}_{0}^{t}+{ \mathcal O }({{\rm{\Sigma }}}_{{pp}}^{2}).\end{eqnarray} \tag{ 69 }$$

Using the zeroth order approximation $p(\tau )\approx {p}_{0}-2F\tau$ in the integral finally yields the first order approximation

$$\begin{eqnarray}&&p(t)\approx {p}_{0}-2{Ft}-\tfrac{g\mathrm{sinh}\mu }{F}({{\rm{\Sigma }}}_{{pp}}(t)\mathrm{sin}({p}_{0}-2{Ft})-{{\rm{\Sigma }}}_{{pp}}(0)\mathrm{sin}({p}_{0})).\end{eqnarray} \tag{ 70 }$$

To approximate ${{\rm{\Sigma }}}_{{pp}}(t)$ in this expression, we substitute the zeroth order approximation for p(t) in the dynamical equation for ${{\rm{\Sigma }}}_{{pp}}$ and integrate to find

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{pp}}\approx \displaystyle \frac{{{\rm{\Sigma }}}_{{pp}}(0)}{1+\displaystyle \frac{g\mathrm{sinh}\mu }{F}(\mathrm{cos}(2{Ft}-{p}_{0})-\mathrm{cos}\;{p}_{0}){{\rm{\Sigma }}}_{{pp}}(0)}.\end{eqnarray} \tag{ 71 }$$

This approximation starts to fail for larger values of the initial momentum width ${{\rm{\Sigma }}}_{{pp}}(0)$, as well as for larger values of μ, as in both cases higher orders would have to be taken into account. These, however, are the cases where the classical approximation as such starts to break down. We demonstrate the onset of this failure in the middle and right panels in figure 3, where we compare the quantum mean value of the momentum to the classical result and the approximation using (70) for the same example as in the left panel, however, with larger initial momentum uncertainty (middle) and a larger value of μ (right panel).

Zoom In Zoom Out Reset image size

Having analysed the behaviour of broad wave packets, let us now turn to the other extreme case where initially only one lattice site is populated. For convenience we choose the central site, i.e., we consider the dynamics of the initial wave function ${c}_{0}(0)=1,{c}_{n\ne 0}(0)=0$. In the Hermitian case this leads to a so-called breathing mode [13], where the centre of the wave packet does not move, while the width oscillates with Bloch frequency. In the non-Hermitian system, this behaviour gets modified, as can be seen in the exact quantum dynamics for an example in the upper left panel in figure 4. we observe that the typical breathing behaviour is strongly biased to one side (depending on the sign of μ), this can be analytically understood from the time-evolution matrix, yielding the time-dependent wave function

$$\begin{eqnarray}&&{c}_{n}(t)={J}_{n}\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{{\rm{i}}(\tfrac{\pi }{2}-{Ft}+{\rm{i}}\mu )n},\end{eqnarray} \tag{ 72 }$$

that is

$$\begin{eqnarray}&&| {c}_{n}(t){| }^{2}={J}_{n}^{2}\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{-2\mu n}.\end{eqnarray} \tag{ 73 }$$

The squared norm of the wave packet can be deduced as

$$\begin{eqnarray}&&P(t)={I}_{0}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right)\end{eqnarray} \tag{ 74 }$$

and the centre is given by

$$\begin{eqnarray}&&\langle \hat{N}\rangle (t)=-\tfrac{2g\mathrm{cosh}\mu }{F}\mathrm{sin}({Ft})\displaystyle \frac{{I}_{1}\left(\tfrac{4g\mathrm{sinh}\mu \mathrm{sin}({Ft})}{F}\right)}{{I}_{0}\left(\tfrac{4g\mathrm{sinh}\mu \mathrm{sin}({Ft})}{F}\right)},\end{eqnarray} \tag{ 75 }$$

where I_ν denotes the modified Bessel functions of first kind.

Zoom In Zoom Out Reset image size

In the Hermitian case it has been shown in [13] that the breathing behaviour can be understood qualitatively in terms of the propagation of a classical ensemble. We shall now demonstrate that this is also true in the non-Hermitian case, and that the classical ensemble dynamics in fact accurately reproduces the dynamics of the centre, the mean momentum, and the squared norm. The initial condition can be interpreted as an extreme case of a Gaussian wave packet, with ${{\rm{\Sigma }}}_{{qq}}=0$ and ${{\rm{\Sigma }}}_{{pp}}\to \infty$. However, we do not expect the classical dynamics to be valid in this case. Thus, we can use the following trick: We can interpret the Fourier transform of our infinitely narrow position Gaussian

$$\begin{eqnarray}&&{\delta }_{n}=\tfrac{1}{2\pi }{\int }_{0}^{2\pi }{{\rm{e}}}^{{\rm{i}}{pn}}{\rm{d}}\;p,\end{eqnarray} \tag{ 76 }$$

as a superposition of plane waves ${{\rm{e}}}^{{\rm{i}}{pn}}$, i.e. Gaussians that are infinitely narrow in space. On the classical side we interpret this superposition as an ensemble of trajectories with initial values ${q}_{0}=0$, ${{\rm{\Sigma }}}_{{pp}}=0$, ${{\rm{\Sigma }}}_{{pq}}=0$, and uniformly distributed p₀. Each of them is then propagated according to equation (61), and weighted with its corresponding squared norm that moves according to equation (62). The resulting propagated classical ensemble is depicted in the upper right panel of figure 4. It can be seen that the breathing feature and the modulation towards the left side are accurately reproduced by the ensemble. Since the classical ensemble does not account for the phases in the quantum superposition features related to interference, such as the horizontal patterns of the propagated beam, are not present in the classical ensemble.

In the bottom panel of figure 4 we show a comparison between the quantum and the classical ensemble averages for the centre and the squared norm of the propagated beam, and observe an excellent agreement. In fact, it can be seen analytically that the ensemble reproduces the exact quantum result for these quantities as follows. The total squared norm of the ensemble is given by the ensemble average of the individual squared norms given in equation (64), that is

$$\begin{eqnarray}\langle P{\rangle }_{{\rm{ensemble}}} & = & \tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }P(t,{p}_{0}){\rm{d}}{p}_{0}\\ & = & \tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }{{\rm{e}}}^{\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\mathrm{sin}({p}_{0}-{Ft})}{\rm{d}}{p}_{0}\\ & = & {I}_{0}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right),\end{eqnarray} \tag{ 77 }$$

where we have assumed that all individual trajectories are initially normalised to one. This result is identical to the quantum squared norm (74) deduced from the time-evolution matrix. The centre of the classical ensemble is similarly given by the ensemble average, weighted by the corresponding squared norms as

$$\begin{eqnarray}\langle q{\rangle }_{{\rm{ensemble}}} & = & \displaystyle \frac{\tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }P(t,{p}_{0})q(t,{p}_{0}){\rm{d}}{p}_{0}}{\langle P{\rangle }_{{\rm{ensemble}}}}\\ & = & -\tfrac{2g\mathrm{cosh}\mu }{F}\mathrm{sin}({Ft})\displaystyle \frac{{I}_{1}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right)}{{I}_{0}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right)},\end{eqnarray} \tag{ 78 }$$

which exactly reproduces the quantum result (75).

We can find the mean value of the quasi momentum from the classical ensemble using a similar argument, as the circular mean of the classical weighted ensemble as

$$\begin{eqnarray}&&\langle p{\rangle }_{{\rm{ensemble}}}=\mathrm{arg}(\langle {{\rm{e}}}^{{\rm{i}}p}{\rangle }_{{\rm{ensemble}}})\end{eqnarray} \tag{ 79 }$$

with

$$\begin{eqnarray}\langle {{\rm{e}}}^{{\rm{i}}p}{\rangle }_{\mathrm{ensemble}} & = & \displaystyle \frac{\tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }P(t,{p}_{0}){{\rm{e}}}^{{\rm{i}}p(t,{p}_{0})}{\rm{d}}{p}_{0}}{\langle P{\rangle }_{\mathrm{ensemble}}}\\ & = & \displaystyle \frac{\tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }{{\rm{e}}}^{\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\mathrm{sin}({p}_{0}-{Ft})}{{\rm{e}}}^{{\rm{i}}({p}_{0}-2{Ft})}{\rm{d}}{p}_{0}}{\langle P{\rangle }_{\mathrm{ensemble}}}\\ & = & \displaystyle \frac{{I}_{1}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right)}{{I}_{0}\left(\tfrac{4g\mathrm{sinh}\mu }{F}\mathrm{sin}({Ft})\right)}\;{{\rm{e}}}^{{\rm{i}}(\tfrac{\pi }{2}-{Ft})}.\end{eqnarray} \tag{ 80 }$$

Thus, we have

$$\begin{eqnarray}&&\langle p{\rangle }_{\mathrm{ensemble}}=\displaystyle \frac{\pi }{2}-{Ft},\end{eqnarray} \tag{ 81 }$$

which also agrees with the exact quantum result that can be obtained using the time-evolution matrix (53).

Let us finally investigate to which extent the ensemble method can be used for intermediately broad wave packets. For this purpose we consider an initial quantum wave packet of the form (65), with $\beta =0.15$, corresponding to ${{\rm{\Sigma }}}_{{pp}}=0.3,{{\rm{\Sigma }}}_{{qq}}=\tfrac{10}{3}$, and vanishing covariance, i.e., ${{\rm{\Sigma }}}_{{pq}}=0$. In figure 5 in the left upper panel we show the renormalised beam propagation for the same parameters as for the previous example, using the exact quantum dynamics. As one might have expected, the observed behaviour is a combination of the dynamics for the two extreme cases of broad and localised wave packets. In the middle panel on the top we show the corresponding classical ensemble dynamics, where the ensemble is distributed according to the Fourier transform of the quantum wave packet, and each trajectory is weighted with its squared norm, as in the previous example. While the classical ensemble fails to correctly describe interference effects as well as the width of the wave packet at the initial state and (correspondingly at times that are multiples of the Bloch period), it describes the dynamical behaviour, and even the width of the beam at most times remarkably well. In comparison we also show the classical propagation of the initial wave packet interpreted as a single classical trajectory with the corresponding initial values of Σ in the right upper panel. Unsurprisingly this classical approximation is not very good, due to the large value of ${{\rm{\Sigma }}}_{{pp}}$. In the bottom panel we show the corresponding comparisons of the dynamics of the centre and the squared norm. We observe that the ensemble and the quantum descriptions agree for both, while the simple classical dynamics gives wrong results for the centre motion, but accurately describes the squared norm.

Zoom In Zoom Out Reset image size

4.2. Purely imaginary coupling constants

As a second example we considered a system with purely imaginary coupling constants ${g}_{1}={g}_{2}={\rm{i}}g$ and a real force $F\in {\mathbb{R}}$. That is, the classical Hamiltonian is given by

$$\begin{eqnarray}&&H=2{\rm{i}}g\mathrm{cos}\;p+2{Fq}.\end{eqnarray} \tag{ 82 }$$

This model was also considered in [9], where an experimental realisation using optical structures has been proposed. The time-evolution matrix for the quantum system can be analytically obtained as [9]

$$\begin{eqnarray}&&{U}_{{nn}^{\prime} }(t)={I}_{n-n^{\prime} }\left(\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{{\rm{i}}(n-n^{\prime} )(\pi -{Ft})-{\rm{i}}2{Ftn}^{\prime} }.\end{eqnarray} \tag{ 83 }$$

Since I_n(x) is real for real arguments and ${I}_{-n}(x)={I}_{n}(x)$, this has the immediate consequence that for an initial state fulfilling ${c}_{-n}={c}_{n}^{* }$ (that is in particular, for an initially real symmetric state), the symmetry is conserved

$$\begin{eqnarray}&&{c}_{-n}(t)={c}_{n}^{* }(t),\end{eqnarray} \tag{ 84 }$$

and therefore we have

$$\begin{eqnarray}&&\langle \hat{N}\rangle (t)=0.\end{eqnarray} \tag{ 85 }$$

This has already been observed for the special case of broad Gaussian wave packets in [9]. We shall now show that this and other features of the quantum dynamics can again be explained on the grounds of the classical dynamics.

The classical equations of motion for the Hamiltonian (82) are given by

$$\begin{eqnarray}&&\dot{p}=-2F-2g\mathrm{sin}\;p\;{{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray}&&\dot{q}=-2g\mathrm{sin}\;p\;{{\rm{\Sigma }}}_{{pq}},\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pp}}=-2g\mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pp}}^{2},\end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{pq}}=-2g\mathrm{cos}\;p\;{{\rm{\Sigma }}}_{{pq}}{{\rm{\Sigma }}}_{{pp}},\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{{qq}}=-2g\mathrm{cos}\;p({{\rm{\Sigma }}}_{{pq}}^{2}-1),\end{eqnarray} \tag{ 90 }$$

$$\begin{eqnarray}&&\dot{P}=-g\mathrm{cos}\;p({{\rm{\Sigma }}}_{{pp}}-4)P.\end{eqnarray} \tag{ 91 }$$

Thus, independently of the width of the wave packet, if the initial covariance of position and momentum vanishes, i.e., ${{\rm{\Sigma }}}_{{pq}}(0)=0$, it stays zero throughout the time evolution. As a consequence also the mean value of the position is constant in time as in the full quantum dynamics. It is interesting to note that this does not imply the conservation of the squared norm, which changes according to the dynamics of the momentum as can be seen for an example in the lower left panel of figure 6.

Zoom In Zoom Out Reset image size

If, however, the covariance is initially non-zero the position performs oscillations in time, induced by the imaginary part of the free dispersion relation. This behaviour is also found in the exact quantum dynamics as demonstrated for a relatively broad initial wave packet with different initial covariances in an example in figure 6, where we observe an excellent agreement between the quantum and the classical dynamics.

For small but non-zero initial momentum uncertainty we again observe periodic modulations of the mean momentum around the linear behaviour predicted by the acceleration theorem, as shown for three examples in figure 7. Using a similar expansion as in the previous example we can approximate these modulations as

$$\begin{eqnarray}&&p(t)\approx {p}_{0}-2{Ft}-\displaystyle \frac{g}{F}{{\rm{\Sigma }}}_{{pp}}(\mathrm{cos}(2{Ft}-{p}_{0})-\mathrm{cos}({p}_{0})),\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{pp}}(t)\approx \displaystyle \frac{{{\rm{\Sigma }}}_{{pp}}(0)}{1-\tfrac{g}{F}(\mathrm{sin}({p}_{0}-2{Ft})-\mathrm{sin}\;{p}_{0}){{\rm{\Sigma }}}_{{pp}}(0)}.\end{eqnarray} \tag{ 93 }$$

This is a good approximation, as long as ${{\rm{\Sigma }}}_{{pp}}(0)$ and g are small, as can be seen in figure 7.

Zoom In Zoom Out Reset image size

Let us finally investigate the ensemble dynamics for an initial state localised in the central lattice site n = 0. The exact quantum evolution is given by

$$\begin{eqnarray}&&{c}_{n}(t)={I}_{n}\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right){{\rm{e}}}^{{\rm{i}}n(\pi -{Ft})},\end{eqnarray} \tag{ 94 }$$

that is, $| {c}_{n}(t){| }^{2}={I}_{n}^{2}\left(-\tfrac{2g}{F}\mathrm{sin}({Ft})\right)$. As discussed above, due to the preserved symmetry, the centre remains stationary at 0, and the wave packet breathes slightly, however, to a much lesser extent than in the previous example or the Hermitian case, as can be seen for an example in the left panel in figure 8. In the middle panel of the same figure we show the classical ensemble propagation, where the dynamics for the individual ensemble trajectories are given by

$$\begin{eqnarray}&&\dot{p}=-2F,\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&\dot{q}=0,\end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&\dot{P}=4g\mathrm{cos}\;p\;P.\end{eqnarray} \tag{ 97 }$$

That is, the momentum follows the acceleration theorem, the position stays constant, and the squared norm evolves as

$$\begin{eqnarray}&&P(t)=\mathrm{exp}\left(\tfrac{4g}{F}\mathrm{sin}({Ft})\mathrm{cos}({Ft}-{p}_{0})\right){P}_{0}.\end{eqnarray} \tag{ 98 }$$

Obviously, since none of the individual trajectories moves, the ensemble prediction is that the whole beam is stationary in the central site, and it fails to describe the breathing of the exact propagation. The squared norm, however, as depicted in the figure on the right is again exactly recovered by the classical ensemble. Similar to the Hatano–Nelson case we find from equation (98)

$$\begin{eqnarray}\langle P{\rangle }_{{\rm{ensemble}}} & = & \tfrac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }{{\rm{e}}}^{\tfrac{4g\mathrm{sin}({Ft})}{F}\mathrm{cos}({Ft}-{p}_{0})}{\rm{d}}{p}_{0}\\ & = & {I}_{0}\left(\tfrac{4g}{F}\mathrm{sin}({Ft})\right),\end{eqnarray} \tag{ 99 }$$

which agrees exactly with the quantum result that can be obtained from the propagator (83). It is interesting to note that in contrast to the Hatano–Nelson model, this always leads to an overall gain. We further find for the circular mean of the quasi momentum

$$\begin{eqnarray}{\langle {{\rm{e}}}^{{\rm{i}}p}\rangle }_{{\rm{ensemble}}} & = & \displaystyle \frac{1}{2\pi {\langle P\rangle }_{{\rm{ensemble}}}}{\displaystyle \int }_{0}^{2\pi }{{\rm{e}}}^{4\tfrac{g}{F}\mathrm{sin}({Ft})\mathrm{cos}({Ft}-{p}_{0})}{{\rm{e}}}^{{\rm{i}}({p}_{0}-2{Ft})}{\rm{d}}{p}_{0}\\ & = & \displaystyle \frac{{I}_{1}\left(\tfrac{4g}{F}\mathrm{sin}({Ft})\right)}{{I}_{0}\left(\tfrac{4g}{F}\mathrm{sin}({Ft})\right)}{{\rm{e}}}^{-{\rm{i}}{Ft}},\end{eqnarray} \tag{ 100 }$$

which also agrees with the exact quantum result, and which yields

$$\begin{eqnarray}&&\langle p{\rangle }_{{\rm{ensemble}}}=-{Ft}.\end{eqnarray} \tag{ 101 }$$

Zoom In Zoom Out Reset image size

We have demonstrated that a quasiclassical description can be used to explain qualitative and even quantitative features of Bloch oscillations in a non-Hermitian single-band tight-binding model. While the examples considered in the present paper are PT-symmetric, the methods developed here can be applied in the same way to more general non-Hermitian systems.

The quasiclassical treatment is a good approximation for wave packets that are narrow in momentum, as the Hamiltonian is anharmonic only in the momentum. For wave packets that are narrow in position space, on the other hand, we have introduced an ensemble method, where we represent the narrow position wave packet as a coherent superpositions of waves that are narrowly localised in momentum space. When evaluating expectation values of observables we approximate this coherent superposition by an incoherent superposition of the mean values of the relevant observable in each of those plane waves. We have observed that this method yields exact results for the norm, the mean position and the mean momentum of wave packets that are initially localised in a single lattice site. In an appendix we provide an analytical derivation of this observation based on the properties of the exact quantum dynamics. The questions of the validity of the results provided by the classical ensemble for other characteristic quantities such as the width of the wave packet, and for initial states that are less localised, as well as to what extent the method is applicable to more general models are interesting topics for future investigations.

We note that the quasiclassical treatment introduced here is a priori only applicable in a single-band approximation, as it is also the case for Hermitian systems. There has been considerable interest in effects of multiband dynamics and non-Hermitian degeneracies between energy bands on Bloch oscillations recently [19, 20]. It is an interesting question to which extend the quasiclassical method could be applied to these cases. Similar to the Hermitian case transitions between the bands would have to be taken into account for example via surface hopping procedures [33].

We have only touched upon the exact quantum dynamics for the two example models considered, but we note that algebraic techniques developed for Hermitian tight-binding systems [26] can be extended to the non-Hermitian case, yielding analytic results for the quantum dynamics, which will be discussed in a forthcoming publication.

E M G acknowledges support from the Royal Society via a University Research Fellowship (Grant. No. UF130339). A R acknowledges support from the Engineering and Physical Sciences Research Council via the Doctoral Research Allocation Grant No. EP/K502856/1.

Here we show that the conjectured interrelation between expectation values of states localised in position and states localised in (quasi)momentum in equation (49) is indeed valid for the present system for operators in the algebra spanned by $\hat{K},{\hat{K}}^{\dagger },\hat{N}$, and the identity, and for wave packets that are initially localised in a single lattice site. Let us denote the squared norm of a state localised initially in the state $| n\rangle$ as well as a state initially in a Bloch state $| \kappa \rangle$ by

$$\begin{eqnarray}&&{P}^{(n)}=\langle \psi | \psi \rangle =\langle n| {\hat{U}}^{\dagger }\hat{U}| n\rangle ,\quad {P}^{(\kappa )}=\langle \psi | \psi \rangle =\langle \kappa | {\hat{U}}^{\dagger }\hat{U}| \kappa \rangle ,\end{eqnarray} \tag{ 102 }$$

respectively, where $\hat{U}=\hat{U}(t)$ is the time evolution operator. The conjectured relation

$$\begin{eqnarray}&&{P}^{(n)}\langle \hat{A}{\rangle }^{(n)}=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{\rm{d}}\kappa \;{P}^{(\kappa )}\langle \hat{A}{\rangle }^{(\kappa )},\end{eqnarray} \tag{ 103 }$$

is equivalent to

$$\begin{eqnarray}&&\langle n| {\hat{U}}^{\dagger }\hat{A}\hat{U}| n\rangle =\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{\rm{d}}\kappa \;\langle \kappa | {\hat{U}}^{\dagger }\hat{A}\hat{U}| \kappa \rangle .\end{eqnarray} \tag{ 104 }$$

Rewriting the left-hand side as

$$\begin{eqnarray}&&\langle n| {\hat{U}}^{\dagger }\hat{A}\hat{U}| n\rangle =\displaystyle \frac{1}{2\pi }{\iint }_{0}^{2\pi }{\rm{d}}\kappa \;{\rm{d}}\kappa ^{\prime} \;{{\rm{e}}}^{{\rm{i}}n(\kappa ^{\prime} -\kappa )}\langle \kappa | {\hat{U}}^{\dagger }\hat{A}\hat{U}| \kappa ^{\prime} \rangle ,\end{eqnarray} \tag{ 105 }$$

immediately shows that equation (104) is valid if ${\hat{U}}^{\dagger }\hat{A}\hat{U}$ is diagonal in the quasimomentum.

This can be conveniently analysed for the present case of a shift algebra formed by the operators $\hat{N}$, $\hat{K}$ and ${\hat{K}}^{\dagger }$ by employing the Wei and Norman exponential product form of the time evolution operator

$$\begin{eqnarray}&&\hat{U}={{\rm{e}}}^{-{\rm{i}}\eta \hat{N}}{{\rm{e}}}^{-{\rm{i}}{\chi }_{1}\hat{K}}{{\rm{e}}}^{-{\rm{i}}{\chi }_{2}{\hat{K}}^{\dagger }},\end{eqnarray} \tag{ 106 }$$

where the coefficients in the exponents are time dependent satisfying the differential equations

$$\begin{eqnarray}&&\dot{\eta }=2F,\quad {\dot{\chi }}_{1}={g}_{1}{{\rm{e}}}^{-{\rm{i}}\eta },\quad {\dot{\chi }}_{2}={g}_{2}{{\rm{e}}}^{{\rm{i}}\eta }\end{eqnarray} \tag{ 107 }$$

with initial conditions $\eta (0)={\chi }_{1}(0)={\chi }_{2}(0)=0$ (see, e.g., [26]). We note that η is real valued for real force constants F, which will be assumed in the following. With ${\hat{U}}^{\dagger }={{\rm{e}}}^{{\rm{i}}{\chi }_{2}^{* }\hat{K}}{{\rm{e}}}^{{\rm{i}}{\chi }_{1}^{* }{\hat{K}}^{\dagger }}{{\rm{e}}}^{{\rm{i}}\eta \hat{N}}$ we find

$$\begin{eqnarray}&&{\hat{U}}^{\dagger }\hat{A}\hat{U}={{\rm{e}}}^{{\rm{i}}{\chi }_{2}^{* }\hat{K}}{{\rm{e}}}^{{\rm{i}}{\chi }_{1}^{* }{\hat{K}}^{\dagger }}{{\rm{e}}}^{{\rm{i}}\eta \hat{N}}\hat{A}{{\rm{e}}}^{-{\rm{i}}\eta \hat{N}}{{\rm{e}}}^{-{\rm{i}}{\chi }_{1}\hat{K}}{{\rm{e}}}^{-{\rm{i}}{\chi }_{2}{\hat{K}}^{\dagger }},\end{eqnarray} \tag{ 108 }$$

which is diagonal in the quasimomentum if ${{\rm{e}}}^{{\rm{i}}\eta \hat{N}}\hat{A}{{\rm{e}}}^{-{\rm{i}}\eta \hat{N}}$ is diagonal. This is obviously the case if the operator $\hat{A}$ is the identity, that is, the norm of a wave packet that is initially localised in space can be exactly deduced from the ensemble method. For $\hat{A}=\hat{K}$ or ${\hat{K}}^{\dagger }$ we can use the operator relation ${{\rm{e}}}^{z\hat{N}}\hat{K}{{\rm{e}}}^{-z\hat{N}}={{\rm{e}}}^{-z}\hat{K}$ or ${{\rm{e}}}^{z\hat{N}}{\hat{K}}^{\dagger }{{\rm{e}}}^{-z\hat{N}}={{\rm{e}}}^{z}{\hat{K}}^{\dagger }$ which shows that (108) is diagonal. For $\hat{A}=\hat{N}$ we have ${{\rm{e}}}^{{\rm{i}}\eta \hat{N}}\hat{N}{{\rm{e}}}^{-{\rm{i}}\eta \hat{N}}=\hat{N}$, which is diagonal in the quasimomentum as already stated in the the discussion of the Hamiltonian in section 2.