A non-Hermitian PT-symmetric kicked top

Non-Hermitian and in particular PT-symmetric quantum systems, that can have either complex conjugate or purely real eigenvalues, have been the subject of extensive investigations over the past two decades. The interest reaches from their spectral theory to technological applications [1, 2]. Yet, apart from a few exceptions [3–6], what has been little investigated hitherto is the interplay between non-Hermiticity and chaos. Since mixed regular-chaotic systems are the norm rather than the exception, it is clearly of great importance to study this interplay to gain deeper insights into the typical behaviour of PT-symmetric quantum systems and their classical counterparts.

Many of the implementations of PT-symmetric systems have used analogous classical wave systems, such as microwaves or classical optics (see e.g. [2] for some examples). Such systems are exceptionally well suited for the investigation of quantum (or more precisely wave) chaos. Recently, however, also quantum mechanical realisations of PT-symmetric systems have been experimentally implemented [7–9]. Here we study the interplay of PT-symmetry and quantum chaos for a toy model based on one of the paradigm models of Hermitian/unitary quantum chaos. While we believe that it is entirely possible to implement this model in a suitable experimental setup, the emphasis of the present study is on the identification of typical behaviour that may be observed in a more readily available experimental system, in the theoretical comfort of a simple model system that is numerically and analytically more accessible.

For this purpose we introduce a PT-symmetric extension of one of the most prominent example systems of traditional quantum chaos: the kicked top. In dependence on the parameters it can show regular, mixed or (almost) entirely chaotic behaviour, and is faithful to random matrix theory [10–12]. More than twenty years after its theoretical inception, it has been experimentally realised using cold atoms [13]. A dissipative version of the kicked top modelling decoherence via Lindblad dynamics has been studied in the literature before [14], and a PT-symmetric version of the non-kicked system has been studied in [15–17]. Here we present the first study of the PT-symmetric kicked system showing a wealth of new behaviour on both the quantum and the classical side, compared to both the non-Hermitian system without PT-symmetry and the Hermitian system. In the classical limit, the dynamics shows a mixture of what may be termed quasi-Hamiltonian chaos and dissipative chaos. This is in fact an example of what has been studied under the name of reversible dynamical systems [18] in classical dynamics.

The paper is organised as follows. We introduce the model of a PT-symmetric kicked top, adding a linear gain-loss profile to the non-kicked part, in section 2, and briefly discuss a possible experimental realisation. In section 3 we derive the classical dynamics using the dynamics of the expectation values of the angular momentum operators, and analyse it in some detail. We provide a very brief introduction to the concept of reversibility in classical maps, which is the classical counterpart of the quantum PT-symmetry. This PT-symmetric/reversible classical dynamics shows stable Poincare dynamics without any sinks or sources, if the gain–loss parameter is not too large. There is no counterpart of this behaviour in non-Hermitian or other open systems without PT-symmetry/reversibility. Furthermore, in our system, we find that the presence of gain and loss leads to an early onset of chaos. Compared to the unitary system, the non-Hermitian PT-symmetric kicked top has an additional classical variable resulting from the norm of the quantum wave function, that can be viewed as a classical intensity. The behaviour of this intensity variable is analysed as a function of the initial state. In section 4 we then study some of the features of the quantum system, starting with the phase-space structure arising from the Husimi representation of the subspaces corresponding to certain lifetimes, which we find to be closely related to the classical intensity as a function of the initial state. Furthermore, the statistics of the nearest-neighbour spacings of the complex quasi energies are analysed to gain an insight into potentially universal features. Recently the spectral statistics of open quantum chaotic systems have attracted renewed interest (see e.g. [19–22]), and new universality classes for non-Hermitian random matrices have been identified in [22]. Here we demonstrate that the nearest-neighbour distributions for our PT-symmetric kicked top in the chaotic regime indeed follow the corresponding universality class of [22]. This universality class is dictated by a symmetry of our model under transposition, that is independent of the PT-symmetry. We show that we can switch to another universality class, by introducing an additional nonlinear kicking term in the Hamiltonian, similar to the switch from GOE to GUE statistics for the standard kicked top analysed in [23]. The resulting behaviour is described by the complex Ginibre ensemble in the large N limit. In the regular regime the (unfolded) nearest-neighbour spacings follow the two-dimensional Poisson distribution. This behaviour is in line with what has been identified over two decades ago by Haake and coworkers for regular and chaotic systems with complex eigenvalues [24–26] and which has recently been reinforced by detailed analysis in [21]. We find that when we exclude eigenvalues too close to the boundary or the symmetry line, the nearest-neighbour spacing statistics is in fact independent of the PT-symmetry, which can be explained due to the latter introducing a global correlation in the spectrum, while the nearest-neighbour correlation is clearly local as has been argued in [21, 22]. We finally consider the eigenvalues statistics of a PT-symmetric generalisation of a triadic Baker map. This system exhibits the same level-spacing characteristics as observed in the PT-symmetric kicked top. We finish with a brief summary in section 5.

We introduce a non-Hermitian kicked top, described by the time-dependent Hamiltonian

$$\begin{equation}\hat{H}\left(t\right)=p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}+\frac{k}{L}{\hat{L}}_{z}^{2}\sum _{n=-\infty }^{\infty }\delta \left(t-n\tau \right),\end{equation} \tag{ 1 }$$

where ${\hat{L}}_{j}$ are angular momentum operators, fulfilling the SU(2) commutation relations

$$\begin{equation}\left[{\hat{L}}_{z},{\hat{L}}_{{\pm}}\right]={\pm}{\hat{L}}_{{\pm}},\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \left[{\hat{L}}_{+},{\hat{L}}_{-}\right]={\hat{L}}_{z},\end{equation} \tag{ 2 }$$

with ${\hat{L}}_{{\pm}}={\hat{L}}_{x}{\pm}\mathrm{i}{\hat{L}}_{y}$. The system has five real-valued parameters p, , γ, k and τ. For γ = 0 we recover the familiar Hermitian kicked top [10–12]. The first two terms describe a linear precession of a quantum angular momentum; the quadratic delta-kicked term leads to the onset of chaos in the classical dynamics. The total angular momentum ${\hat{L}}^{2}={\hat{L}}_{x}^{2}+{\hat{L}}_{y}^{2}+{\hat{L}}_{z}^{2}$ is a conserved quantity. Thus, one can restrict the analysis to a (2L + 1)-dimensional sector of the Hilbert space belonging to the eigenvalue L(L + 1) of ${\hat{L}}^{2}$. In the classical picture this restricts the system to a spherical phase space. The kicked top has been experimentally realised using the spin angular momentum of a single ¹³³Cs atom as reported in [13]. Another interpretation of the system is an ensemble of bosonic particles populating two modes (either in an external potential or using internal degrees of freedom) where the interaction between the particles is periodically switched on for infinitesimally short times [27], which could be experimentally realised in cold atoms similar to the experiments reported in [28].

For non-vanishing γ the Hamiltonian is no longer Hermitian, that is, the dynamics is not unitary, and the norm is not conserved. The γ term describes a linear gain/loss profile along the L_z direction. In the bosonic two-mode context this could be realised (up to an imaginary energy shift) by a loss of particles from one of the modes. A related model without explicit time-dependence, where the quadratic part is constantly switched on, has been studied theoretically in detail in [16, 17]. Experimentally atom losses from lattices can be controlled with high precision [29, 30]. They are naturally described by Lindblad dynamics, and an effective non-Hermitian dynamics can be achieved by post selection (see, e.g. [31]), as has been the case in the recent experimental realisations of quantum PT-symmetric systems [9, 32]. In this spirit, the kicked top (1) could be realised with atoms in a double well trap, with induced losses from one of the sites, and a time-dependent interaction strength, or equivalently, time-dependent barrier height and loss rate [27].

There are various definitions of PT-symmetry in use throughout the literature, that differ in small details. Here we use the definition of [33]: a non-Hermitian Hamiltonian is called PT-symmetric if it satisfies [PT, H] = 0 where T is an antilinear operator and P is a linear involution that commutes with T. T is usually set to be the time-reversal operator, which, in our case fulfils T² = +1 (sometimes referred to as even PT-symmetry). A PT-symmetric Hamiltonian Ĥ with even PT-symmetry thus satisfies Ĥ = PĤ*P, where the asterisk denotes complex conjugation. As a consequence Ĥ has either real eigenvalues when the corresponding eigenstates are also eigenstates of the operator PT, or otherwise the eigenvalues come in complex conjugate pairs. For = 0 for any given time the Hamiltonian (1) is PT-symmetric, with

$$\begin{equation}PT:\quad {\hat{L}}_{x}{\mapsto}{\hat{L}}_{x},\quad {\hat{L}}_{y}{\mapsto}{\hat{L}}_{y},\quad {\hat{L}}_{z}{\mapsto}-{\hat{L}}_{z},\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \mathrm{i}\to -\mathrm{i}.\end{equation} \tag{ 3 }$$

The investigation of kicked systems is most easily facilitated by the study of the Floquet operator, describing the time-evolution of a state over one time period, $\hat{F}=\hat{U}\left({t}_{0},{t}_{0}+\tau \right)$. Here, we choose a time-interval symmetric around the kick, that is, we consider a free evolution with ${\hat{H}}_{0}=p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}$ for time length $\frac{\tau }{2}$, followed by an instantaneous kick with ${\hat{H}}_{1}=\frac{k}{L}{\hat{L}}_{z}^{2}$, and then another free evolution for time $\frac{\tau }{2}$. That is, the Floquet operator is given by

$$\begin{equation}\hat{F}={\mathrm{e}}^{-\mathrm{i}\left(p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}\right)\frac{\tau }{2}}\enspace {\mathrm{e}}^{-\mathrm{i}\frac{k}{L}{\hat{L}}_{z}^{2}}\enspace {\mathrm{e}}^{-\mathrm{i}\left(p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}\right)\frac{\tau }{2}}.\end{equation} \tag{ 4 }$$

The choice of the time interval ensures that the Floquet operator is indeed PT-symmetric with respect to the PT operator above for = 0. Here we follow [34] and define a quantum map as PT-symmetric if it fulfills the relation $\hat{F}=P{\left({\hat{F}}^{{\ast}}\right)}^{-1}P$. Denoting the eigenvalues of a PT-symmetric quantum map by λ = re^−iθ, the eigenvalues are either unimodular, when the corresponding eigenstates are also eigenstates of the PT operator, or come in pairs ${\lambda }_{1}={\left({\lambda }_{2}^{{\ast}}\right)}^{-1}$, i.e., ${r}_{1}=\frac{1}{{r}_{2}}$ and θ₁ = θ₂, in which case the PT operator acting on one of the eigenvectors yields the other. That is, the quasienergies E = i ln(λ) are either real or come in complex conjugate pairs, just as the eigenvalues of time-independent PT-symmetric Hamiltonians. Note that for another choice of the time interval the Floquet operator would be isospectral to the one considered here and still be PT-symmetric, however, with a non-trivial P operator.

In what follows we shall study the occurrence of chaos in the classical dynamics and its fingerprints in the quantum system for the PT-symmetric quantum map (4). For this purpose, we shall first introduce the classical map, arising in the limit of large Hilbert space dimension L → ∞, and analyse the resulting dynamics, before we return to its quantum features.

To derive the classical dynamics we consider the quantum evolution of the expectation values of the three angular momentum operators in the limit L → ∞. We will deduce the mapping of the expectation values of the angular momentum operators over one period, by expressing

$$\begin{equation}{\langle {\hat{L}}_{j}\rangle }_{n+1}=\frac{\langle {\psi }^{\left(n+1\right)}\vert {\hat{L}}_{j}\vert {\psi }^{\left(n+1\right)}\rangle }{\langle {\psi }^{\left(n+1\right)}\vert {\psi }^{\left(n+1\right)}\rangle }\end{equation} \tag{ 5 }$$

as functions of the ${\langle {\hat{L}}_{j}\rangle }_{n}$, where

$$\begin{equation}\vert {\psi }^{\left(n+1\right)}\rangle =\hat{F}\vert {\psi }^{\left(n\right)}\rangle ={\hat{F}}^{\left(n+1\right)}\vert {\psi }^{\left(0\right)}\rangle ,\end{equation} \tag{ 6 }$$

is the quantum state after n + 1 evolution periods.

Let us first consider the free evolution with the Hamiltonian ${\hat{H}}_{0}=p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}$. The equations of motion for the expectation values of the angular momentum operators can for example be deduced from the 2 × 2 representation of the system, which, due to the linearity in the ${\hat{L}}_{j}$ provides the full dynamics for all values of the total angular momentum. The time evolution of the expectation values of the rescaled angular momentum expectation values ${s}_{j}=\langle {\hat{L}}_{j}\rangle /L$ for the free evolution is given by [16, 17]

$$\begin{align}\hfill \begin{aligned}\hfill {\dot {s}}_{x}& =-{\epsilon}{s}_{y}-\gamma \enspace {s}_{z}{s}_{x},\hfill \\ \hfill {\dot {s}}_{y}& ={\epsilon}{s}_{x}-p{s}_{z}-\gamma \enspace {s}_{z}{s}_{y},\hfill \\ \hfill {\dot {s}}_{z}& =p{s}_{y}+\gamma \enspace \left(1-{s}_{z}^{2}\right).\hfill \end{aligned}\end{align} \tag{ 7 }$$

In the Hermitian case γ = 0 this describes a rigid rotation around the axis ${\left(p,\enspace 0,\enspace {\epsilon}\right)}^{T}$.

Some examples of the phase-space trajectories for non-vanishing γ are depicted in figure 1. For more details on the dynamics see, e.g., [17]. Figure 1 depicts both the dynamics on the sphere (in the top row) as well as the projection onto a flat phase space in (canonical) cylindrical coordinates s_z and $\phi =\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{{s}_{y}}{{s}_{x}}\right)$ (in the bottom row). In the PT-symmetric case, γ ≠ 0, = 0, there are two distinct phases of motion for subcritical and supercritical values of γ, respectively, where the critical value is given by γ_crit = p. For γ < p, an example of which is shown in the left panel of figure 1, the motion is a combination of a rotation around the x-axis and a Lorentz boost around the z-axis, with two stable elliptic fixed points that, for increasing values of γ move from the opposite poles at s_x = ±1 along the equator, approaching s_x = 0, where they meet at the critical value. Note that while there are no sinks and sources in this evolution the dynamics is not area preserving, but has balanced amounts of contraction and expansion over the phase sphere. For supercritical increasing values of γ > γ_crit the fixed points separate and travel along the great circle with s_x = 0 towards the north and south pole. The southern fixed point is a source and the northern fixed point is a sink of the dynamics. An example of this motion is depicted in the middle panel of figure 1. This classical behaviour is directly related to the PT-symmetry and its breaking in the quantum system [17]. Without the PT-symmetry, for ≠ 0 the elliptic fixed points present for the Hermitian system turn into a sink and a source for arbitrarily small values of γ, and the rigid rotation around the ${\left(p,\enspace 0,\enspace {\epsilon}\right)}^{T}$-axis turns into a spiralling phase-space flow from source to sink, where again the sink and the source move closer to the poles of the sphere for increasing values of γ. An example of the resulting trajectories is depicted in the right panel of figure 1.

Zoom In Zoom Out Reset image size

To deduce the phase-space map generated by the free motion over half a period τ/2 we need to integrate the dynamical equation (7). For the PT-symmetric case = 0, that explicitly yields

$$\begin{equation}{\tilde {s}}_{x}=\enspace \enspace \frac{1}{{\Gamma}}{s}_{x}\end{equation} \tag{ 8 }$$

$$\begin{equation}{\tilde {s}}_{y}=\enspace \enspace \frac{1}{{\Gamma}}\left(\left(\frac{{p}^{2}}{{\omega }^{2}}\enspace \mathrm{cos}\left(\frac{\omega \tau }{2}\right)-\frac{{\gamma }^{2}}{{\omega }^{2}}\right){s}_{y}-\frac{p}{\omega }\enspace \mathrm{sin}\left(\frac{\omega \tau }{2}\right){s}_{z}+\frac{p\gamma }{{\omega }^{2}}\enspace \mathrm{cos}\left(\frac{\omega \tau }{2}\right)-\frac{p\gamma }{{\omega }^{2}}\right)\end{equation} \tag{ 9 }$$

$$\begin{equation}{\tilde {s}}_{z}=\enspace \enspace \frac{1}{{\Gamma}}\left(\enspace \frac{p}{\omega }\enspace \mathrm{sin}\left(\frac{\omega \tau }{2}\right){s}_{y}+\mathrm{cos}\left(\frac{\omega \tau }{2}\right){s}_{z}+\frac{\gamma }{\omega }\enspace \mathrm{sin}\left(\frac{\omega \tau }{2}\right)\right),\end{equation} \tag{ 10 }$$

with

$$\begin{equation}{\Gamma}=\frac{{p}^{2}}{{\omega }^{2}}-\frac{{\gamma }^{2}}{{\omega }^{2}}\enspace \mathrm{cos}\left(\frac{\omega \tau }{2}\right)+\frac{p\gamma }{{\omega }^{2}}\left(1-\mathrm{cos}\left(\frac{\omega \tau }{2}\right)\right){s}_{y}+\frac{\gamma }{\omega }\enspace \mathrm{sin}\left(\frac{\omega \tau }{2}\right){s}_{z},\end{equation} \tag{ 11 }$$

and $\omega =\sqrt{{p}^{2}-{\gamma }^{2}}$. It is straightforward to verify that this indeed reduces to a rigid rotation by an angle $\frac{p\tau }{2}$ around the x-axis, for γ = 0.

The full kicked top map is composed of the mapping resulting from the free motion followed by the mapping arising from the kick, which is not altered by the non-Hermiticity of the free evolution, and then another free motion of half a period. The kicked part of the Floquet operator alone yields the map [11]

$$\begin{align}\hfill {\tilde {s}}_{x}& =\mathrm{cos}\left(2k{s}_{z}\right){s}_{x}-\mathrm{sin}\left(2k{s}_{z}\right){s}_{y}\hfill \\ \hfill {\tilde {s}}_{y}& =\mathrm{sin}\left(2k{s}_{z}\right){s}_{x}+\mathrm{cos}\left(2k{s}_{z}\right){s}_{y}\hfill \\ \hfill {\tilde {s}}_{z}& ={s}_{z},\hfill \end{align} \tag{ 12 }$$

that is, a torsion around the z-axis. In what follows we shall set the period τ to unity.

The PT symmetry of the quantum system results in the classical map F = F_freeF_kickF_free being reversible [18]. A classical map F map is called reversible if there exists an involution G (i.e. applying G twice yields the identity) such that

$$\begin{equation}F{\circ}G{\circ}F=G\end{equation} \tag{ 13 }$$

which is equivalent to

$$\begin{equation}{F}^{-1}=G{\circ}F{\circ}G,\end{equation} \tag{ 14 }$$

i.e. under the transformation of the phase space by G the map (that is the time evolution) is inversed. Not surprisingly, the involution G for our system is the classical counterpart of the P operator, that is the reversing of the z-direction G = P: s_z ↦ −s_z.

There are many consequences of reversibility, regarding symmetry lines, the occurrence of pseudo conservative features such as periodic orbits and KAM tori, as well as dissipative features such as attractors, and universal scaling behaviours [18]. An exhaustive analysis of the classical dynamics goes well beyond of the scope of the present paper, however, and will be reserved for another study. Here we shall content ourselves with gaining a first overview of the dynamical features of our classical map.

In the Hermitian case γ = 0 an increasing kicking strength k leads to the onset of chaos and for large kicking strength k ≫ p the dynamics is almost entirely chaotic. Here we shall investigate how the non-Hermiticity changes this behaviour, in particular in the PT-symmetric case = 0.

Note that even for infinitesimal values of γ the dynamics is no longer area preserving. Since the kicking part of the dynamics remains Hamiltonian, the rate of contraction and expansion does not depend on the value of the kicking strength. Some examples of Poincare plots of the phase-space dynamics for = 0, p = 2, and relatively small γ and k are depicted in figure 2. We observe that the non-Hermiticity seems to enhance the chaos, but does not lead to the appearance of sinks and sources in the dynamics for the small values depicted here. This behaviour is due to the PT-symmetry/reversibility of the system, and an inclusion of a small symmetry breaking parameter would lead to a fully contracting phase-space flow with a single point sink.

Zoom In Zoom Out Reset image size

The enhancement of chaos in the classical phase-space dynamics observed in the Poincare sections is confirmed by the analysis of the Lyapunov exponent of the dynamics, depicted in false colours for p = 2, k = 1 and three different values of γ in figure 3. A positive Lyapunov exponent indicates exponential sensitivity to small variations in the initial conditions, that is, chaos. In the examples depicted, while most of the phase space is associated to a vanishing Lyapunov exponent, corresponding to regular motion, for vanishing γ, the region of positive Lyapunov exponent grows significantly with increasing γ and for γ = 0.5 already a large proportion of the phase space is covered by a chaotic sea.

Zoom In Zoom Out Reset image size

The most prominent fixed point in the Poincare sections in figure 2 is the one that originates from the elliptic fixed point at ϕ = 0 for vanishing γ and k. It moves along the equator with increasing γ, much as it would have done for k = 0. In fact, its position is unaffected by the value of k, provided k is not so large that the fixed point has bifurcated, which can be seen as follows. It is a consequence of the reversibility of the system that any possible elliptic fixed point has to lie on the line s_z = 0. In addition, at an elliptic fixed point the linearisation of the map must be area preserving, i.e. the Jacobian determinant must be one. Now that Jacobian determinant, as mentioned above, is independent of the (conservative) kicking and thus the parameter k. Therefore, the only possible elliptic fixed points of the system are the two fixed points of the free evolution, that are situated at s_z = 0 and ${\phi }_{1}=-\mathrm{arctan}\left(\frac{\gamma }{p}\right)$ and ${\phi }_{2}=-\pi -\mathrm{arctan}\left(\frac{\gamma }{p}\right)$, respectively. However, with increasing k the fixed point on the left side in our phase space plots bifurcates into a two-cycle and a hyperbolic fixed point in a chaotic band. For not too large k, as γ increases further, the bifurcation is reversed, and the stable fixed point on the left reappears, as can be seen for example in figure 4. The regions around the two fixed points are separated from each other by increasing numbers of layers of chains of small regular islands and chaotic bands, in the typical KAM fashion.

Zoom In Zoom Out Reset image size

For larger, but still subcritical values of γ < p the prominent central fixed point moves further along the equator and finally bifurcates into a sink and a source for a small volume surrounding it. An example of the Poincare section just above this bifurcation, and a magnification of the area around the central fixed point, are depicted in figure 4 in the left and middle panels. At this value of γ also other of the countless fixed points of the system have turned into sinks for tiny phase-space volumes, as is revealed by a zoom into the Poincare plot, depicted in the right panel of the same figure. To generate these figures both the forwards and backwards time-evolution have been plotted, to get a more complete picture of the phase-space flow.

For k = 1, for further increasing γ the basin of attraction of the sink in the central region grows. At another critical value of γ, finally this region of attraction reaches the chaotic sea, and now the volume of the latter is also attracted by the sink, which by now is the only sink of the flow. An example of the resulting Poincare map and a magnification of the area around the sink and source are depicted in the top left and middle panels of figure 5. The right panel on the top of the same figure shows the same magnified area as the figure in the middle, however, for the backwards propagation, revealing the source corresponding to the sink. The sink attracts all of the phase-space volume but the small volume in the remaining regular island at s_z = 0 and $\phi =-\pi -\mathrm{arctan}\left(\frac{\gamma }{p}\right)$. This can be seen in the bottom panels of the same figure, that depict the distribution of 20 000 phase-space points that were initially chosen from a uniform distribution on phase space and propagated for 100, 500, and 10, 000 iterations, respectively. It can be clearly seen that while the points that started out in the regular island remain within this volume, all other points are eventually mapped into the sink.

Zoom In Zoom Out Reset image size

For even larger values of γ the remaining regular island shrinks further and moves further along the equator approaching the hyperbolic fixed point between the sink and the source that moves towards it along the equator from the other side. At the critical value γ = p the stable island finally vanishes when it meets the hyperbolic fixed point, and the whole phase-space volume is attracted to the single remaining sink. While the details of the dynamics in dependence on γ depend strongly on k, independently of the value of k there are no conservative features in the dynamics for γ > p. This is related to the fact that the dynamics does no longer allow for elliptic fixed points above this critical value. Of course, depending on k the transition away from pseudo conservative phase space portraits can happen for much smaller γ already.

To provide a better idea of typical phase-space portraits we depict a number of further examples in figure 6 for two different values of γ and increasing values of k. Note that for intermediate values of k as in the middle column, most of the phase-space volume is attracted to the sink on the northern hemisphere, much as in figure 5. For larger values of k we clearly see a spiral shaped fractal structure emerging. Depending on the exact values of γ and k this structure can be transient while for long times the phase-space volume is attracted to a smaller sink, or it can be a fully developed strange attractor.

Zoom In Zoom Out Reset image size

In the Hermitian case γ = 0 increasing values of the kicking strength k lead to enlarged chaotic seas spreading over the phase sphere that eventually almost cover the sphere for sufficiently large values of k. An example of the Poincare section for p = 2 and k = 7 is depicted on the left in figure 7. Introducing gain and loss, γ ≠ 0, in this case leads to the appearance of spiral-shaped attractors in the phase space, as can be seen in the middle and right panels of the same figure. These spirals appear to be strange attractors. That is, over time the dynamics continues to move around the spiral in a random manner, and zooms into the Poincare plots, as depicted for example in figure 8, indeed reveal a fractal structure of the attractor.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shape of these strange attractors can be easily understood as being structured around the skeleton of the orbits arising from the single period map, that is, the free motion distorted by the torsion of the kicking dynamics and then followed by another free motion. The left panels in figure 9 depict the image of the great circle under this transformation on the sphere as well as on the projection for the same values of p and k as in figure 7. We have chosen γ = 0 for this purpose, since the deformation of the spiral for small values of γ is negligible. This demonstrates clearly, that the seed of the strange attractor actually lies within the Hermitian system, but is made visible only by the non-Hermiticity. Using the image of the great circle under this transformation to describe the strange attractor, we can easily deduce that the winding number of the spiral is given by $\frac{4k}{\pi }$, since the north pole is rotated by 2kπ in the positive direction, and the south pole is rotated by the same angle in the negative direction by the kicking torsion. We can further deduce that the 'centre' of the spiral, is located at ϕ = −π (i.e. s_x = 0), and ${s}_{z}=\mathrm{cos}\left(\frac{p}{2}\right)$. In the middle and right panels of figure 9 this dependence of the winding number and orientation of the attractor is demonstrated using two different sets of parameters. The top row shows the Poincare sections of the map, while the bottom row shows the image of the great circle under the transformation of the kicking torsion followed by the free motion.

Zoom In Zoom Out Reset image size

Note that in contrast to the behaviour for smaller k, where an arbitrarily small value of the PT-symmetry breaking parameter leads to appearance of a single sink in the phase-space dynamics, these strange attractors are stable and also appear in the non-Hermitian system without PT-symmetry. We have verified for a large range of parameters that the motion indeed continues to erratically move around the attractor rather than going into a sink for large time scales with or without PT-symmetry.

Finally, the classical limit of non-Hermitian systems has an additional variable on top of the phase-space coordinates, that describes the behaviour of the norm of the quantum state in the semiclassical limit,

$$\begin{equation}{n}_{j+1}{:=}\langle {\psi }^{\left(j+1\right)}\vert {\psi }^{\left(j+1\right)}\rangle =\langle {\psi }^{\left(j\right)}\vert {\hat{F}}^{{\dagger}}\hat{F}\vert {\psi }^{\left(j\right)}\rangle .\end{equation} \tag{ 15 }$$

We have already encountered the map for this variable in the 2 × 2 representation as

$$\begin{equation}{n}_{j+1}={\Gamma}\left({s}_{x},{s}_{y},{s}_{z}\right){n}_{j},\end{equation} \tag{ 16 }$$

with Γ defined in equation (11). This additional parameter can be viewed as an intensity of the classical phase-space trajectory, which growth and decays locally exponentially under the free dynamics proportionally to the momentary z-component, but is not affected by the (unitary) kick. In general the intensity is a complicated function of the iteration time for most initial states. To provide an overview of its behaviour over time, we show false-colour plots of the final intensity as a function of the initial phase-space variable for different parameter values and different numbers of iterations in figure 10. We observe that as the number of iterations is increased finer and finer structures emerge. There is an obvious correspondence to the phase-space portraits. In particular in the bottom panel, belonging to the case of the strange attractor, the fractal structure is very clearly visible in the final intensity after increasing numbers of iterations. Note however, that the structure visible in the norm dynamics does not correspond to the attractor, but the repeller of the motion, that would be obtained from a backwards evolution in time.

Zoom In Zoom Out Reset image size

We further find that for large values of k for which the unitary map is chaotic, even for arbitrarily small values of γ the fractal repeller is clearly visible in the iterated intensity, while it is not visible in the dynamics, with the bare eye. This can be seen for example in the middle panel of figure 13 where the classical intensity after three iterations as a function of the initial state is depicted for p = 2, k = 10, and γ = 0.01.

After having uncovered some of the rich behaviour of the classical dynamics of this PT-symmetric chaotic system, let us now turn our attention to the quantum system.

Let us now consider how the classical features observed in the previous section manifest themselves in the full quantum description. In particular, we will investigate the spectral behaviour of the Floquet operator and the phase-space structures associated to its eigenspaces. The quantum system is very sensitive to non-Hermitian perturbations, and an accurate numerical diagonalisation is problematic for medium sized values of L if the gain–loss parameter is too large. This sensitivity to perturbations is related to the presence of high-order exceptional points in the non-kicked system [15], and is indeed typical for non-Hermitian and non-unitary (more precisely non-normal) operators. To illustrate the issue we show two examples of selected contours of the resolvent norm (where we take the two-norm, computed using the singular value decomposition) of the quantum Floquet operator for a relatively small matrix size of 41 (L = 20) in figure 11.

Zoom In Zoom Out Reset image size

The eigenvalues of an operator are defined as the singularities of the resolvent norm, that is, points in the complex plane at which the resolvent norm diverges, i.e., values λ at which ${\Vert}\hat{A}-\lambda \mathbb{I}{\Vert}=0$. The theory of pseudospectra [35] deals with the concept of values λ which are almost eigenvalues in the sense that ${\Vert}\hat{A}-{\lambda }_{{\epsilon}}\mathbb{I}{\Vert}={\epsilon}$, or equivalently ${\Vert}{\left(\hat{A}-{\lambda }_{{\epsilon}}\mathbb{I}\right)}^{-1}{\Vert}=\frac{1}{{\epsilon}}$ for arbitrarily small , that is, points at which the resolvent norm is large but finite. For normal operators the resolvent norm is very sharply peaked around eigenvalues and the points for which the resolvent norm is larger than $\frac{1}{{\epsilon}}$ only differ from the actual eigenvalues on the order of , which confirms the intuition that a pseudo-eigenvalue λ is indeed close to an eigenvalue [35]. Yet, this intuition turns out to be misleading for operators that are not normal, such as our Floquet operator. In fact, our system, due to the presence of many exceptional points at which eigenvectors coalesce, is a particular notorious example of this phenomenon. This can be seen in the right plot in figure 11 where the level curves at which the resolvent norm equals 10⁸ to 10¹¹ are plotted for a relatively small matrix size of 41 × 41, and moderate values of p = 2, k = 1, and γ = 1. The level curves extend over a large range of quasi energy values and we can see that solutions to the equation ${\Vert}\hat{F}-{\lambda }_{{\epsilon}}\mathbb{I}{\Vert}=1{0}^{-9}$ can differ from each other, and thus the true eigenvalue on the order of 10. This phenomenon is hugely amplified for larger matrix sizes and quickly leads to perturbations on the order of machine accuracy leading to pseudo eigenvalues that are largely uncorrelated with the true eigenvalues. In fact, the right conclusion from this behaviour may be to follow Trefethen's advice [35] and not to consider eigenvalues at all for our system, as they may say very little about the systems behaviour. However, for the present purpose we stick with the analysis of the spectral behaviour of our quantum system to identify typical features of quantum–classical correspondence, but we are forced to confine ourselves to small values of γ for which the spectral problem is better behaved, as illustrated for the same values of L, p, and k, as before, but γ = 0.1 in the left panel of figure 11.

Let us begin with a small value of k for which the classical system is in the regular regime. The left panel in the top row of figure 12 shows the quasi energies in the complex plane for p = 2, and k = 0.5 and a very small loss–gain parameter γ = 0.01 for L = 200. The corresponding regular classical dynamics is visualised by a Poincare section in the middle panel of the top row. While the dynamics does not show any sinks and sources, the intensity is not conserved for most initial conditions, but varies slowly, due to the small loss–gain parameter, which introduces an asymmetry in the trajectories. This is reflected in the quantum quasi energies, a handful of which lie visibly off the real axis. The PT-symmetry of the system is clearly obvious from the symmetry of the quasi energies with respect to the real axis. Since the system is only mildly non-Hermitian most eigenstates are still approximately orthogonal to each other, and their phase-space distributions are organised according to the structure of the classical phase-space dynamics. This is illustrated in the remaining four panels in figure 12, which show the Husimi distributions of four different Floquet states.

Zoom In Zoom Out Reset image size

In the regime which is chaotic in the unitary case, on the other hand, a small non-Hermitian perturbation has a very different effect. Figure 13 illustrates the case p = 2, k = 10 and γ = 0.01, again for L = 200. In this case almost all the quasi energies (depicted in the complex plane in the left panel of the figure) are complex. Further, the Floquet matrix is far from normal, and most eigenvectors strongly overlap with other eigenvectors. Nevertheless, one can see some structures reminiscent of the classical attractor and repeller in the Husimi functions of individual Floquet states, three of which are depicted in the middle and right panel on the top and the left panel in the bottom row. To obtain a better insight into the quantum phase-space properties, rather than analysing individual eigenstates, it is useful to consider the Husimi representation of certain invariant subspaces, the Husimi–Schur representation as introduced in [36]. In the middle panel in the bottom row of figure 13 we depict the Husimi representation of the invariant subspace of states with growing norm, that is, we use the Schur vectors of the Floquet operator where the diagonal elements are sorted according to the imaginary part of the quasi energies, and average over the Husimi distributions of the Schur vectors belonging to the positive imaginary parts. The resulting picture shows a striking resemblance to the false colour plot arising from the classical intensity as a function of the initial phase-space position, which is depicted after 3 iterations in the right panel. This correspondence takes the place of the familiar quantum phase-space localisation on classical Poincare structures in phase space for the unitary case for strongly open quantum systems. Similar observations have been made for other types of openings in chaotic systems, for example in [36–38]. It is an interesting open question how the eigenspaces of the quantum system are organised in phase space in regions where the real energy structures and the norm dynamics are competing and neither dominates. This question, however goes beyond the scope of the present paper.

Zoom In Zoom Out Reset image size

This classical phase-space structure is also visible in the quantum dynamics, as demonstrated in figure 14, which depicts the Husimi representation of an initial eigenstate of ${\hat{L}}_{z}$ after different numbers of iterations. The spiral structure is clearly visible in the iterated states. After a large number of iterations all that remains is the fastest growing state. Indeed for the present example after approximately 20 iterations the state has approached the final configuration. This can be seen from the similarity of the Husimi representation of the state after 30 iterations, depicted in the lower right panel in figure 14 and the Husimi function of the fastest growing state, depicted in the lower left panel in figure 13.

Zoom In Zoom Out Reset image size

Let us finally turn to the spectral statistics of our PT-symmetric kicked top. In the Hermitian case the level spacing statistics of the quasi energies of the kicked top follow GOE statistics in the chaotic regime, and show a Poisson distribution in the regular regime, in line with the Bohigas–Giannoni–Schmit conjecture. In the non-Hermitian case, the quasi-energies become complex, and a generalisation of level spacings needs to be considered. Here we follow the approach in [24–26] and analyse the distributions of Euclidean nearest neighbour distances in the complex plane, after unfolding the spectra. For the unfolding we make use of the fact that the quasi energy spectrum is approximately uniform in the real parts, and use the integrated staircase function of the imaginary parts to approximate the underlying smooth distribution and use this to unfold the imaginary parts, while leaving the real parts unscaled. To obtain a suitable amount of quasienergies we combine the data from a number of realisations with slightly altered parameters, similar to the common practice in the Hermitian case. In figure 15 we depict the nearest-neighbour distance distributions in the regular regime, for k varying in 0.05 increments between k = 0.1 and k = 1, for L = 200 and a small γ = 0.01, where we disregarded quasi energies close to the boundary to avoid boundary effects. The nearest-neighbour distance distribution in the chaotic case, for values of k = 9 to k = 10 (in increments of 0.02) and γ = 0.01 is depicted in the top row of figure 16. For comparison we also plot the two-dimensional Poisson distribution curve (green) that is expected for independent random numbers in the plane, and associated to regular behaviour for complex eigenvalues [24–26]. This distribution is identical to the GOE distribution

$$\begin{equation}P\left(s\right)=\frac{\pi }{2}s{\mathrm{e}}^{-\frac{\pi }{4}{s}^{2}}.\end{equation} \tag{ 17 }$$

The second curve (magenta), shows the nearest neighbour distribution for eigenvalues of Gaussian random matrices with the constraint  = Â^T, which we numerically obtained from the diagonalisation of 10 000 matrices of the same size as for the kicked top (401 × 401). This has been identified as one of three universality classes for non-Hermitian Hamiltonians in [22], which are speculated to replace the three universality classes of Hermitian systems given by Dyson's Gaussian ensembles. Indeed, our Hamiltonian and Floquet operator are invariant under matrix transposition and we observe an excellent agreement of the data in the regular regime with the two-dimensional Poisson curve, and with the Gaussian transpose invariant random matrix model in the chaotic regime. The more familiar universal nearest-neighbour distance distribution class for non-Hermitian matrices that has been identified in the spectra of various physical systems in the past [21, 24–26] is depicted by the cyan curves in figure 16. The nearest-neighbour spacing distribution is given by that of the eigenvalues of the Ginibre ensembles in the large N limit as [26]

$$\begin{equation}P\left(s\right)=C\enspace {\mathrm{lim}}_{N\to \infty }\left(\prod _{n=1}^{N-1}{e}_{n}\left({\left(Cs\right)}^{2}{\mathrm{e}}^{-{\left(Cs\right)}^{2}}\right)\right)\sum _{n=1}^{N-1}\frac{2{\left(Cs\right)}^{2n+1}}{n!{e}_{n}\left({\left(Cs\right)}^{2}\right)},\end{equation} \tag{ 18 }$$

with

$$\begin{equation}{e}_{n}\left(x\right)=\sum _{m=0}^{n}\frac{{x}^{m}}{m!},\end{equation} \tag{ 19 }$$

and the numerical constant (arising from demanding normalisation and that ⟨s⟩ = 1)

$$\begin{equation}C\approx 1.1429\cdots \enspace .\end{equation} \tag{ 20 }$$

This behaviour can be recovered if we introduce an extra kicking term proportional to ${\hat{L}}_{y}^{2}$ in our model to yield

$$\begin{equation}\hat{F}={\mathrm{e}}^{-\mathrm{i}\left(p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}\right)\frac{\tau }{2}}\enspace {\mathrm{e}}^{-\mathrm{i}\frac{k}{L}{\hat{L}}_{z}^{2}}{\mathrm{e}}^{-\mathrm{i}\frac{\eta }{L}{\hat{L}}_{y}^{2}}\enspace {\mathrm{e}}^{-\mathrm{i}\left(p{\hat{L}}_{x}+\left({\epsilon}+\mathrm{i}\gamma \right){\hat{L}}_{z}\right)\frac{\tau }{2}},\end{equation} \tag{ 21 }$$

much as in the unitary case, where this extra factor changes the level spacing behaviour from GOE to GUE [23]. The level-spacing distribution for the same values of L, p, k, and γ and with η = 1 are depicted in the bottom row of figure 16, where we observe an excellent agreement between the data and the prediction (18) from the Ginibre ensemble.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have further analysed the spectra of a PT-symmetric generalisation of another quantum chaotic model showing GOE statistics in the Hermitian case, a triadic Baker map. There are multiple quantisations of the Baker map, and many investigations of open versions (see e.g. [39–41]). Here we use a standard quantisation in which the unitary triadic Baker map is given by the Floquet operator (see, e.g. [39] and references therein)

$$\begin{equation}\hat{\mathcal{B}}={F}_{N}^{-1}\left(\begin{matrix}\hfill {F}_{N/3}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {F}_{N/3}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {F}_{N/3}\hfill \end{matrix}\right),\end{equation} \tag{ 22 }$$

where F_N is a discrete Fourier transform,

$$\begin{equation}{\left({F}_{N}\right)}_{mn}=\frac{1}{\sqrt{N}}\enspace {\mathrm{e}}^{-\frac{2\pi \mathrm{i}}{N}\left(m+\frac{1}{2}\right)\left(n+\frac{1}{2}\right)}; m,n\in \left\{0,1,\dots ,N-1\right\}.\end{equation} \tag{ 23 }$$

This map is parameter free, the classical counterpart is fully chaotic, and the quantum quasi energies display random matrix statistics. We generalise this to a non-unitary PT-symmetric map (with ${\mathcal{T}}^{2}=\mathbb{I}$) by introducing a parameter γ ∈ (0, 1) that amplifies phase-space trajectories in one third of the phase space, leaves the central third unchanged, and introduces a loss (of equal strength to the gain) in the remaining third, yielding the map

$$\begin{equation}{\hat{\mathcal{B}}}_{PT}={F}_{N}^{-1}\mathrm{diag}\left(\gamma {F}_{N/3},{F}_{N/3},\frac{1}{\gamma }{F}_{N/3}\right).\end{equation} \tag{ 24 }$$

Figure 17 shows the statistics of the nearest-neighbour distances of the quasi energies for the PT-symmetric Baker map with matrix size 3000 averaged over different values of γ ranging in 0.01 steps from γ = 0.5 to γ = 0.55. The Baker map also falls in the class of transpose invariant matrices since the similarity transformation ${F}_{N}^{\frac{1}{2}}{\hat{\mathcal{B}}}_{PT}{F}_{N}^{-\frac{1}{2}}$ indeed brings it to a transpose-invariant form, and we observe a close agreement with the corresponding random matrix predictions.

Zoom In Zoom Out Reset image size

Preliminary studies of a PT-symmetric generalisation of an M-mode Bose–Hubbard model with intermediate interaction strength and unit filling factor also show similar behaviour. This system is known to show the typical quantum classical correspondence of chaotic systems in its many-body and mean-field correspondence (see, e.g. [42–44]). This offers an interesting opportunity to observe features of dissipative quantum chaos in many-body quantum systems, as losses in cold atoms in optical lattices are an active area of research and can be well controlled experimentally (see, e.g. [29, 45, 46]).

Thus, we conclude that the statistics of the nearest-neighbour spacings in the complex plane for a number of PT-symmetric models are indeed universal. This particular statistics, however in the bulk of the spectrum is independent of the PT-symmetry of the system. It would be interesting to investigate whether other spectral statistics differ between PT-symmetric and more general non-Hermitian models. A natural starting point for this would be the investigation of the statistical features of the real eigenvalues alone, that do not have a counterpart in generic non-Hermitian systems. This, however, would require a much larger data set, and goes beyond the scope of the present paper.

In summary, we have introduced a simple PT-symmetric generalisation of a quantum kicked top. We have derived the classical analogue of the dynamics and found that it leads to chaotic phase-space dynamics without sinks or sources for a wide parameter range, and a mixture of stable orbits and dissipative behaviour for larger gain-loss parameters, and strange attractors similar to conventional dissipative chaotic systems when the corresponding unitary system is deep in the chaotic regime. Similar behaviour has been discussed before in the theory of classically reversible dynamical systems, which indeed are the natural counterpart of PT-symmetric quantum systems. The classical dynamics arising here has an additional degree of freedom, which can be interpreted as an intensity, that shows intricate dynamical behaviour. In the full quantum system, the nature of the classical counterpart is reflected in the eigenvalue behaviour as well as in quantum phase-space structures. Finally, we have considered the spectral statistics of the quasi energies, and found universal behaviour. We have verified the same behaviour for a PT-symmetric generalisation of a triadic Baker map. Much remains to be explored. To give just two examples, a detailed characterisation of the rich features of the classical dynamics as well as the investigation of possible universal spectral features that are due to the PT-symmetry are interesting tasks for future studies.

The authors thank Martin Kratky for useful input into the classical dynamics, and Joseph Hall and Simon Malzard for kindly providing a code for the calculation of the Lyapunov exponent, which has been adapted for the present purpose. Further thanks goes to the anonymous referees who have helped to improve this study considerably with their exceptionally insightful questions. EMG acknowledges support from the Royal Society (Grant no UF130339) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement no 758453). SMN acknowledges support from an EPSRC DTA grant (Grant no EP/M506345/1).