Chaos in circuit QED: decoherence, localization and nonclassicality

While chaos is well defined in classical mechanics, the degree to which it can manifest itself in quantum mechanics remains a problem of fundamental interest [1]. For example, classically chaotic systems exhibit exponential sensitivity to initial conditions, something that is absent from quantum dynamics which is explicitly unitary. Moreover, the discrete spectrum of closed quantum systems implies quasi-periodicity which contradicts the irregular behaviour seen in classical chaotic systems. Theoretical work has indicated that the inclusion of decoherence, on the other hand, can cause a quantum system to become exponentially sensitive to initial state fluctuations [2].

On the experimental front, the advent of the laser cooling of atoms has allowed the study of a number of quantum chaotic phenomena in simple and highly controllable settings. Of particular relevance to this paper we highlight the observation of dynamical localization of momentum in a kicked quantum rotor [3] made from atoms in a frequency modulated optical lattice, where classical (chaotic) diffusion is suppressed by quantum interference, and chaos assisted tunnelling [4, 5], which was observed using atoms in an amplitude modulated lattice. In the latter system the tunnelling is dynamical in nature in the sense that it occurs in phase space between islands of regular motion, allowing the system to pass through the classically insurmountable dynamical 'barriers' provided by the boundaries between integrable and chaotic regions. The kicked quantum top has also been realized using atomic spins [6]: using continuous weak measurement techniques, the full quantum dynamics in phase space can be reconstructed as a function of time. Some of these experimental setups have also been used to investigate the effects of decoherence [7, 8], and, as expected, it was found that in each case the system reverted to classical behaviour.

Driven by the promise of technological applications, continuing progress in the manufacture of man made and hybrid quantum systems also provides significant opportunities for exploring basic quantum dynamics under controllable conditions [9]. Within the realm of this progress the idea of quantum simulators has become reality [10–12]. While quantum simulation of closed systems has been the main focus thus far, the extension to open systems has been considered [13]. In our view these new quantum systems naturally lend themselves to the study of chaotic quantum dynamics subject to decoherence.

In this work we explore quantum chaos in a circuit quantum electrodynamics (QED) setup comprised of a superconducting qubit resonantly coupled to a single mode of a high-Q transmission line resonator [14]. The presence of the two-level qubit means that this is a manifestly quantum system in the sense that it does not possess a simple classical limit. The parameters of the corresponding Hamiltonian can be externally controlled to a large extent, and even the strength of the decoherence can be manipulated by various means. In order to achieve chaotic dynamics we consider the case where the resonator is externally pumped by a drive oscillator with adjustable amplitude and phase. By changing either the qubit–field coupling [15], the pump amplitude or the pump modulation frequency, the system can be varied between being regular or chaotic. In particular, in the ultrastrong coupling regime [15], and at vigorous enough pumping, the system becomes chaotic in the sense that the semi-classical equations of motion corresponding to a classical electromagnetic field (but retaining the quantum nature of the qubit), become chaotic. Furthermore, in the Born–Oppenheimer approximation (BOA) where the excited state of the qubit is adiabatically eliminated, we find that the driven circuit QED system mimics a quantum Duffing oscillator [16]. The Duffing oscillator is an example of a well-studied system whose classical limit is known to be chaotic.

In the chaotic regime we shall see that an initially localized phase space distribution of the electromagnetic field rapidly broadens and the Wigner distribution builds up seemingly irregular sub-Planck structures characteristic of chaos. In comparison to how a localized phase space distribution evolves under the corresponding semi-classical equations of motion, the broadening in the fully quantum evolution is much more dramatic leading to a breakdown of the quantum-classical correspondence. One possible route to recover the classical dynamics is to try to localize the state by weakly measuring the position quadrature of the field. This induces decoherence and non-unitary time evolution [17]. One of the main results of our work concerns the effect of such decoherence on the present chaotic bipartite system. In general, it has been shown that for quantum chaotic systems that possess a clear classical limit, decoherence is a possible candidate to explain how chaos, in a classical sense, emerges in quantum mechanics. In particular, decoherence can prevent rapid spreading of an initially localized state and the quantum-classical correspondence remains valid for long times [2]. However, when the system does not support a classical limit, like here, the situation becomes more complex [18].

In the present model the effects of quantum fluctuations in the two-level system are substantial and we cannot ascribe it a direct classical counterpart. As a consequence, we will demonstrate that the semi-classical predictions are not recovered by introducing measurement induced decoherence of the oscillator subsystem (electromagnetic field) in our model. We attempt to understand this lack of correspondence by showing that fundamentally non-classical characteristics, such as qubit–field entanglement, survive the decoherence. In the optical regime, the effects of a weak coupling to an environment or a measurement device is conveniently captured in terms of a master equation [19, 20]. In an unravelling approach to solving the master equation, the Lindblad terms constitute imaginary state decay together with stochastic quantum jumps. Averaging over all these quantum trajectories results in the evolved state of the system. Single stochastic trajectories, which can be interpreted as single experimental measurement realizations, are shown to support localization. As will be explained further below, it is important to distinguish this localization from dynamical localization [1, 3, 21, 22] which is a quantum coherent effect that is destroyed by decoherence not generated by it [7]. The localization phenomenon we discuss here only appears in the states of single stochastic trajectories and not in the full quantum probability distribution (i.e. the average over an ensemble of measurements). We argue that signatures of chaotic behaviour in circuit QED can be seen in the photon statistics of the transmission line output field.

This paper is organized as follows. The next section introduces the model system and the decoherence stemming from the position measurement. In section 3 we derive a set of semi-classical equations of motion where the field variable is treated at a mean-field level while a quantum treatment is maintained for the qubit. The results of the full quantum model are presented in section 4. We especially focus on dynamics in phase space which can be visualized using Wigner and Husimi distributions, including localization appearing due to 'selective measurements', and methods to detect signatures of quantum chaos. Finally, we conclude in section 5. In the appendix we present pictures of typical examples of the residuals of the phase space probability distributions in order to give a better idea of how the individual x and p probability distributions look under the combined effects of unitary evolution and measurement.

Our physical model consists of a superconducting qubit coupled to a driven transmission line resonator as, for example, in the experiments described in [14, 23]. The undriven setup is described by the Rabi Hamiltonian [24] (ℏ = 1)

$$\begin{equation} \hat{H}_R=\omega \left(\frac{\hat{p}^2}{2}+\frac{\hat{x}^2}{2}\right)+\frac{\Omega }{2}\hat{\sigma }_z+\sqrt{2}g\hat{x}\hat{\sigma }_x, \end{equation} \tag{ 1 }$$

where ω and Ω are the mode and qubit transition frequencies respectively, g the effective qubit–light interaction strength, $\hat{x}$ and $\hat{p}$ the two field quadratures of the resonator mode, and $\hat{\sigma }_x=|2\rangle \langle 1|+|1\rangle \langle 2|$ and $\hat{\sigma }_z=|2\rangle \langle 2|-|1\rangle \langle 1|$ are the Pauli matrices acting on the two qubit states |1〉 and |2〉. We will work with dimensionless variables where ℏω sets the energy scale, i.e. from now on ω = 1. The above Hamiltonian can be microscopically derived from a minimal coupling model where the two-level, single mode, and dipole approximations are imposed, and the 'self-energy' of the field is neglected, see [25]. In the rotating wave approximation (RWA) one recovers the celebrated Jaynes–Cummings model [24, 26]. In traditional microwave and optical cavity QED setups, g ≪ 1 such that this approximation is well justified [27], while recent advances in circuit QED make it possible to reach regimes where the RWA breaks down [28]. When one goes beyond the RWA limit, especially for the ultrastrong coupling regime $g>\sqrt{\Omega }/2$, the photon population $\bar{n}_0$ of the system ground state increases rapidly ($\bar{n}_0\sim g^2$). For the traditional configuration, namely a two-level system interacting with a single resonator mode, this build-up of photon population $\bar{n}_0$ in the ground state is mainly due to the fact that the electromagnetic self-energy in the Hamiltonian (1) has been left out [29]. This problem can be circumvented by following the timely suggestion made in [15]: by adding an external drive for the qubit one obtains an effective model of the Rabi form, but where the qubit–field coupling is scaled with the pump amplitude while the self-energy is unaffected, and, therefore, the two terms can be tuned independently. Hence, the self-energy can be neglected even in the ultrastrong coupling regime of the effective Rabi model. This fact is important for us here because we find that chaotic structures appear most definitively in the ultrastrong coupling regime.

Making use of the Z₂ parity symmetry of the Rabi Hamiltonian, the model was recently shown to be quasi-solvable over its entire parameter range [30]. This integrability is somewhat surprising since the model does not possess a continuous symmetry. Moreover, integrability typically implies absence of chaos and the dynamics is therefore regular for all energies and parameters [31]. However, as we will demonstrate, when the system is externally driven it may become chaotic for certain time-dependent drive amplitudes. The driven Rabi Hamiltonian is given by

$$\begin{eqnarray} &&\hat{H}_{dR} = \displaystyle {\hat{a}^\dagger \hat{a}+\frac{\Omega }{2}\hat{\sigma }_z+g (\hat{a}^\dagger +\hat{a})\hat{\sigma }_x} +\eta (t)(\hat{a}^\dagger {\rm e}^{-{\rm i}\omega _d t}+\hat{a}{\rm e}^{{\rm i}\omega _dt}), \nonumber \\ \end{eqnarray} \tag{ 2 }$$

where we have introduced the photon annihilation and creation operators $\hat{a}=(\hat{x}+{\rm i}\hat{p})/\sqrt{2}$ and $\hat{a}^\dagger =(\hat{x}-{\rm i}\hat{p})/\sqrt{2}$, respectively. The drive amplitude will be taken to be

$$\begin{equation} \eta (t)=\eta _0\cos (\omega _ct), \end{equation} \tag{ 3 }$$

i.e. the amplitude is modulated in time at angular frequency ω_c between the maximum/minimum values of ±η₀. The amplitude modulation frequency should not be confused with the central frequency ω_d of the drive field upon which it imposes side bands. Returning to the quadrature representation, we have

$$\begin{eqnarray} &&\hat{H}_{dR} = {\frac{1}{2}(\hat{p}+\sqrt{2}\eta (t)\sin (\omega _dt))^2+\frac{\hat{x}^2}{2}+\frac{\Omega }{2}\hat{\sigma }_z} \nonumber \\ &&{(\sqrt{2}g\hat{\sigma }_x+\sqrt{2}\eta (t)\cos (\omega _dt))\hat{x}-\frac{\eta ^2(t)\sin ^2(\omega _dt)}{2}}, \end{eqnarray} \tag{ 4 }$$

which will be our starting Hamiltonian.

To gain further insight, it is convenient to perform an adiabatic diagonalization where the Hamiltonian is expressed in the eigenvalue basis of the two-level matrix $V(x)=\Omega \hat{\sigma }_z/2+gx\hat{\sigma }_x$ with corresponding adiabatic potentials

$$\begin{equation} V_{{\rm ad}}^\pm (x,t)=\frac{x^2}{2}+\sqrt{2}\eta (t)\cos (\omega _dt)x\pm \sqrt{\frac{\Omega ^2}{4}+2g^2x^2}. \end{equation} \tag{ 5 }$$

The BOA consists of neglecting any coupling between the resulting adiabatic eigenstates [32]. If we let η₀ = 0 we note that the lower potential $V_{{\rm ad}}^\pm (x)$ has a double-well structure whenever $g>\sqrt{\Omega }/2$. The driving induces a back and forth rocking of the double-well. This type of driven double-well system mimics the Duffing oscillator [16]. However, for non-zero ω_d the driving also induces a 'shaking' of the momentum. From equation (4) it follows that this shaking effect can be accorded a time-dependent gauge potential $\hat{A}(t)=\sqrt{2}\eta (t)\sin (\omega _dt)$ giving rise to a synthetic electric field $\hat{E}=-\partial _t\hat{A}(t)$. While interesting, in this paper we shall not discuss this gauge aspect further. As will be described in the next section, the phenomena we are interested in stem from the rocking of the potential and this synthetic electric field will not change our arguments.

Finally we want to include the impact of quadrature measurements [33]. In the Schrödinger picture and within the Born–Markov approximation a standard approach for describing decoherence is given by considering a master equation where the irreversible processes are attributed to a set of Lindblad operators [20]. Instead of considering field and qubit decay, we consider decoherence arising from the measurement of the $\hat{x}$-quadrature of the field. Thus, we work with a Markovian master equation of the form

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{\rho }={\rm i}[\hat{\rho },\hat{H}_{dR}]+\mathcal {L}[\hat{x}]\hat{\rho }, \end{equation} \tag{ 6 }$$

where $\hat{\rho }$ is the density operator for the full system, and the Lindblad operator

$$\begin{equation} \mathcal {L}[\hat{x}]\hat{\rho }=\kappa (2\hat{x}\hat{\rho }\hat{x}^\dagger -\hat{x}^2\hat{\rho }-\hat{\rho }\hat{x}^2)=-\kappa [\hat{x},[\hat{x},\hat{\rho }]] \end{equation} \tag{ 7 }$$

accounts for the irreversible loss of coherence due to the measurement of the quadrature $\hat{x}$ [33]. The information loss rate is given by the parameter κ. Intuitively, a large κ is expected to produce a localization effect upon the quantum state. However, if the measurement back-action becomes too large it will inevitably affect the variable conjugate to the measured one. Thus, in order to achieve a localization effect, κ should be in between the two extremes. In the present bipartite model this picture will, however, be very different. Note that the form of $\mathcal {L}[\hat{x}]$ implies only loss of coherence and no dissipation. This will allow us to identify any localization effect in phase space as being due to decoherence rather than a reduction of the available phase space due to energy loss. We point out that our results are not unique to decoherence of the form of (7) but other decoherence channels like $\mathcal {L}[\hat{a}^\dagger \hat{a}]$ give similar behaviour (as verified numerically). For pure photon decay $\mathcal {L}[\hat{a}]$ which also induces dissipation, the quantum-classical crossover in a Jaynes–Cummings model with a Kerr medium was studied in [34] in terms of collapse-revivals.

A further remark before proceeding: whenever the system at hand is the subject of decoherence, information about it is irreversibly lost leading to an effective nonlinear theory that can become chaotic in a classical sense, i.e. for short times nearby initial states move apart exponentially with a rate set by the maximum Lyapunov exponent [35]. However, in an analysis like ours, any irreversible process in quantum mechanics derives from the effective nature of the theory. Normally, it is assumed that the measurement device is macroscopic so that information about the system is distributed over a large number of degrees-of-freedom and a reversed flow of information back into the system is very unlikely. This justifies the Markovian assumption of infinitely fast decay of correlations within the measurement device [9, 20]. It is nevertheless an approximation and by waiting long enough information will flow back into the system (quantum recurrences). This waiting time will typically be very long in comparison to any experimentally relevant time scale so that the physics will appear to be nonlinear in the experiment.

At large pumping, η₀ > 1, we replace the mode variables by their corresponding c-numbers as obtained in a coherent state ansatz. Then, the resulting effective Hamiltonian acts on the qubit sub-space alone, while its parameters depend on the boson field and must be determined self-consistently. In such a mean-field approach we neglect any quantum correlations between the qubit and the field, i.e. the qubit remains in a pure state and its presence does not imply decoherence of the large amplitude resonator field. This allows us to derive a closed set of equations of motion, and thereby characterize the underlying semi-classical dynamics, i.e. whether it is regular, chaotic, or mixed.

Since the qubit is decoupled from the field, we can describe it via two variables. After denoting the qubit state as $|\chi \rangle =[\beta _1{\rm e}^{-{\rm i}\phi _1}\,,\,\beta _2{\rm e}^{-{\rm i}\phi _2}]^{\rm T}$, where β_j and ϕ_j (j = 1, 2) are real, we introduce the (conjugate variables) inversion $Z=\beta _1^2-\beta _2^2$ and relative phase Δ_ϕ = ϕ₁ − ϕ₂. The semi-classical equations of motion become

$$\begin{eqnarray} &&{\dot{x}= p+\sqrt{2}\eta (t)\sin (\omega _dt),} \nonumber \\ &&{\dot{p}=- x-g\sqrt{2}\sqrt{1\!-Z^2}\cos (\Delta _\phi )-\sqrt{2}\eta (t)\cos (\omega _dt),} \nonumber \\ &&\dot{Z}=g\sqrt{2}x\sqrt{1-Z^2}\sin (\Delta _\phi ), \nonumber \\ &&{\dot{\Delta }_\phi =\frac{\Omega }{2}-g\sqrt{2}x\cos (\Delta _\phi )\frac{Z}{\sqrt{1-Z^2}}}, \end{eqnarray} \tag{ 8 }$$

where the dot denotes a time-derivative. We have left out the decoherence induced Langevin forces [19] in this mean-field analysis. As a result, the semi-classical dynamics is independent of κ which reflects the fact that the quadrature measurement only induces decoherence and no dissipation. From equations (8), we find the semi-classical Hamiltonian

$$\begin{eqnarray} &&H_{\rm scl}(x,p,Z,\Delta _\phi ) = \displaystyle {\frac{p^2}{2}+\frac{x^2}{2}+p\sqrt{2}\eta (t)\sin (\omega _dt)} \nonumber \\ &&+\,x\sqrt{2}\eta (t)\cos (\omega _dt)\\ &&+\,gx\sqrt{2}\sqrt{1-Z^2}\cos (\Delta _\phi )+\frac{\Omega }{2}Z. \nonumber \end{eqnarray} \tag{ 9 }$$

We stress that the mean-field approximation is only applied to the boson degree-of-freedom, in equation (8) the spin is treated quantum mechanically.

From equation (8) we directly identify a trivial stable fixed point given by x_ss = p_ss = 0 whenever η₀ = 0 (i.e. undriven system). At g² = Ω/4, there is a Pitchfork bifurcation where the fixed point turns unstable and two new stable fixed points arise (g² > Ω/4), $x_{ss}=\pm g\sqrt{2-\frac{\Omega ^2}{8g^4}}$. These correspond to the minima of the double-well potential $V_{{\rm ad}}^-(x)$. In an adiabatic picture, the driving induces an effective time-dependent 'detuning' between the two double-well minima. For moderate driving, this dynamics can be treated as repeated Landau–Zener transitions [36]. For large pumping η₀ (typically |η₀| > |g|), however, the double-well structure may be lost at the extreme of the driving, i.e. when |η(t)| attains its maximum. In this case, the lower adiabatic potential oscillates between having one or two local minima.

In this paper we use the semi-classical equations of motion to perform an analysis of the dynamics within the truncated Wigner approach (TWA) [37]. That is, we simulate N_cl ≫ 1 trajectories of the above equations (8) with initial conditions for x and p randomly picked according to a Gaussian probability distribution corresponding to a coherent state |α〉 with amplitude $\alpha =x_{ss}/\sqrt{2}$. The qubit initial state is taken Δ_ϕ(t = 0) = π and Z(t = 0) = 0. This implies that the initial state predominantly populates the right minimum of the lower adiabatic potential $V_{{\rm ad}}^-(x)$. The set of differential equations (8) are solved using the fourth order Runge–Kutta algorithm as modified by Gear's method, which is especially suited to stiff equations.

Throughout this work we will consider the resonant situation Ω = ω_d = ω = 1, the pump amplitude η₀ = 3, amplitude modulation frequency ω_c = 1.5, and effective atom field coupling g = 1.5. This choice means that we are well into the ultrastrong coupling regime $g>\sqrt{\Omega }/2$, and, as we shall see below, η₀ = 3 is large enough to guarantee classical chaos. It should be pointed out that the results are not sensitive to the particular parameter values as long as the dynamics is chaotic. In fact, by fine tuning the drive parameters, chaotic evolution is also attainable in the 'normal phase' $g<\sqrt{\Omega }/2$ where the lower adiabatic potential has a global minimum at x = 0.

With ω_c = 3ω_d/2, the Hamiltonian is periodic with period T = 4π. We consider stroboscopic maps of the phase space distribution of the resonator mode at times t = T and 4T. As can be seen in figure 1, after one period the distribution still shows some visible structures, while after four periods almost all such structures are washed out. For this choice of initial state and system parameters, the dynamics is completely chaotic lacking any islands of regular motion [35].

Zoom In Zoom Out Reset image size

For a time-independent system, chaotic dynamics typically smears out the distribution over the accessible phase space energy shells. We note that the corresponding effect in our time-dependent system occurs at rather short time scales, e.g. the classical period of the harmonic oscillator T_cl = 2π in scaled dimensionless units. It is clear that the stroboscopic maps of figure 1 show different characteristics compared to those of the quantum Duffing oscillator [38]. This results from the additional degree-of-freedom of the spin. The effective phase space is here four dimensional instead of two dimensional as it is for the Duffing oscillator and therefore the typical stretching and squeezing structure of the distribution [38] is not visible in the stroboscopic maps shown in figure 1 since they are projections onto the two dimensional xp-plane.

The master equation (6) gives the evolution of the quantum state $\hat{\rho }$ under the influence of a certain form of decoherence which can be interpreted as a continuous position measurement. Extracting $\langle \hat{x}\rangle$ implies measuring the field quadrature for an ensemble of equally prepared initial states [27]. The statistical average of such homodyne measurement sequences gives $\langle \hat{x}\rangle$. Similarly, when solving (6) numerically it is practical to transform the problem into one of solving a sequence of independent stochastic Schrödinger equations, each one of them resulting in a solution |ψ_i(t)〉 (note that this is a pure state). In this unravelling method any expectation value becomes $\langle \hat{A}\rangle =\mathrm{Tr}\, [\hat{A}\hat{\rho }(t)]=\lim _{N\rightarrow \infty }\frac{1}{N}\sum _i^N\langle \psi _i(t)|\hat{A}|\psi _i(t)\rangle$. Each $\langle \hat{A}\rangle _i=\langle \psi _i(t)|\hat{A}|\psi _i(t)\rangle$ can be viewed as one of the measurement sequences and every photon detection (a 'click' in the detector) is represented by a single 'quantum jump' in the stochastic Schrödinger equation. It is important to appreciate the distinction between the outcome from a single measurement sequence and one deriving from an ensemble of measurements. We note that both are of experimental relevance.

For the numerics we employed the split-operator method [39] combined with the quantum Monte Carlo method [40]. The split-operator method relies on factorizing, for short times δt, the evolution operator into a real-space and a momentum part. The corresponding two propagators are applied in real and momentum space respectively. The method becomes exact in the limit δt → 0. In general, the inherent time scales set the constraints on the size of δt. The final time is t = Kδt, where K is the number of steps used in the propagation. We varied δt in order to check the convergence of our extracted results. The split-operator method allows us to simulate the coherent evolution, while for decoherence we apply the quantum Monte Carlo method, in which N different trajectories |ψ_i(t)〉 are simulated and where for each simulation i we periodically interrupt the coherent propagation and randomly choose either the absence or occurrence of a quantum jump characterizing a photon detection [27]. This interruption is done K times, i.e. it agrees with the time-steps used for the coherent propagation, and we have verified that the jump probability P_jump ≪ 1 in each step. Note that a quantum jump in the present homodyne scheme is represented by multiplying the state |ψ(t)〉 by $\hat{x}$ (and not by $\hat{a}$ as in a normal detection scheme where the photon intensity leaking out of the cavity is measured). As an initial state we will consider a coherent state predominantly occupying one of the two minima of $V_{{\rm ad}}^-(x,t=0)$. Such a state can be approximated as |ψ_±(t = 0)〉 = | ± α〉|Θ_±〉, where $\alpha =x_{ss}/\sqrt{2}$ is a real coherent state amplitude and $|\Theta _\pm \rangle =[\pm 1\,,\,1]^{\rm T}/\sqrt{2}$ is the initial qubit state. The ± sign represents the left or right potential minimum, respectively. This state approximates the ones used for the initialization of the TWA considered in the previous section. Notice that the initial state represents a disentangled field–qubit state.

A Markovian master equation of the Lindblad form may imply (depending on the actual form of the Lindblad terms) dissipation and decoherence. In a driven system, dissipation may balance the pumping such that a steady state containing a finite average number of photons is approached [19]. The main result of decoherence, on the other hand, is to suppress the off-diagonal terms of the density operator $\hat{\rho }$ [41]. This diagonalization can hinder the evolution, similar to the quantum Zeno effect. For chaotic systems, the role of reservoir induced decoherence has been discussed in numerous works, see for example [2, 38, 42, 43]. Decoherence stemming from quantum measurements has been explored as well, see [17, 33]. In the case of the quantum Duffing oscillator, it has been verified that the expectation values $\langle \hat{A}\rangle _{t}$ obtained from the corresponding classical model agreed much better with those of the quantum model when decoherence is included than those given by the closed quantum system [38]. Moreover, the expectation values can become sensitive to small fluctuations in the initial states when decoherence is included, a typical characteristic for chaotic dynamics. At the same time, it has long been known that decoherence is the main source for loss of 'quantumness', e.g. quantum entanglement or negativity of the Wigner function [9, 20]. It might therefore be expected that decoherence helps in obtaining a proper quantum-classical correspondence for chaotic systems. While the issue of decoherence in chaotic systems has been rather well studied for cases in which the system has an obvious classical limit, for systems with no clear classical limit little is known [18]. The situation is even more complex in bi- or multipartite systems where entanglement may persist between the different subsystems even in the presence of decoherence [44].

4.1. Phase space analysis

Quantum phase space distributions are a practical tool for achieving a deeper intuition for quantum dynamical systems [19, 45]. The idea is to visualize the quantum evolution in a phase space picture. The Wigner distribution for a quantum state $\hat{\rho }(t)$ (in our case the resonator boson mode) is defined as [45]

$$\begin{equation} W(x,p,t)=\frac{1}{\pi }\int {\rm d}y\,\langle x-y/2|\hat{\rho }(t)|x+y/2\rangle \, {\rm e}^{{\rm i}py}. \end{equation} \tag{ 10 }$$

Because it does not have to be positive everywhere in phase space, W(x, p, t) is not a proper probability distribution. Typically, negativity of W(x, p, t) is identified with non-classical phenomena such as quantum coherent superpositions. The Wigner distribution is normed to unity, ∫∫dx dp W(x, p, t) = N₋(t) + N₊(t) = 1, where N_±(t) measures the negative (−) and positive fraction (+) of the distribution. Note that normalization does not imply |N_±(t)| < 1. The marginal distributions of W(x, p, t) give the real and momentum space probability distributions, i.e. they are positive definite. Another characteristic of W(x, p, t) is that it allows for sub-Planck structures [46] which is another indication that the Wigner function cannot be regarded as a formal probability distribution. As a counterpart to classical phase space distributions, the Husimi Q-function seems more of a natural choice. The Q-function is positive definite and in contrast to the Glauber P-functions, it does not become singular. It is defined as [19, 45]

$$\begin{equation} Q(\tilde{x},\tilde{p},t)=\frac{1}{\pi }\langle \alpha |\hat{\rho }(t)|\alpha \rangle , \end{equation} \tag{ 11 }$$

where |α〉 is a coherent state with amplitude $\alpha =\alpha _r+{\rm i}\alpha _i\equiv \sqrt{2}(\tilde{x}+{\rm i}\tilde{p})$ (α_r and α_i are hence both real). It follows that the Husimi Q-function is a Gaussian averaged Wigner function [47],

$$\begin{equation} Q(\tilde{x},\tilde{p},t)=\frac{1}{\pi }\int\!\!\!\int {\rm d}x \, {\rm d}p\,{\rm e}^{-(x-\tilde{x})^2-(p-\tilde{p})^2}W(x,p,t). \end{equation} \tag{ 12 }$$

This smoothening not only extinguishes the negativity, but also prevents sub-Planck structures from forming [48]. Interestingly, despite this averaging, because of the overcompleteness of coherent states the Q-function still contains all the information there is to know about the state $\hat{\rho }(t)$. While the Q-function has the appealing property of being positive definite, its marginal distributions do not reproduce the correct real and momentum space probability distributions. Note that the phase space distributions are defined for general states $\hat{\rho }(t)$, with no restriction on purity. This applies especially to us where we are interested in the distributions for the field. The density operator for the field $\hat{\rho }_f(t)$ is typically in a mixed state even without coupling to the environment since it is obtained from the full density operator by a trace over the qubit's degrees-of-freedom, $\hat{\rho }_f(t)=\mathrm{Tr}_\mathrm{qu}[\hat{\rho }(t)]$.

For a chaotic system like the present one, we expect the phase space distributions to show an ergodic collapse in which they spread out over the accessible phase space. Even at resonant driving, ω_d = Ω = 1, the model appears to have a relatively bounded energy shell due to the oscillating pump amplitude η(t) and the coupling to the qubit. After the collapse, a seemingly irregular phase space structure builds up. In a closed quantum system, the above behaviour is typical for quantum thermalization [49]. How the phase space structures develop in a chaotic open system is, however, not well explored. In our model we cannot talk about quantum thermalization since the system is both driven and exposed to decoherence.

We have numerically calculated W(x, p, t) and $Q(\tilde{x},\tilde{p},t)$ for different decay rates κ. The results are shown in figures 2 and 3 respectively, while figure 4 gives a zoom-in of a single Planck cell. Comparing the figures, it is apparent that the sub-Plank structures in the Wigner distribution (ℏ = 1) do not occur in the Q-function. Furthermore, as expected, decoherence enters as a smoothening of the Wigner distribution, and for large κ the sub-Planck structures are completely lost in the Wigner distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Instead of deriving the distributions directly from the state $\hat{\rho }_f(t)$, they can alternatively be obtained from the corresponding Fokker–Planck equations [19]. This approach gives us additional insight into the phase space dynamics. The full bipartite Q-function [19], $Q(z,\alpha )=\frac{2}{\pi ^2(1+|z|^2)}\langle z,\alpha |\hat{\rho }(t)|z,\alpha \rangle$, where |z〉 is a spin coherent state with complex amplitude z [50], obeys the equation

$$\begin{equation} \partial _tQ(z,\alpha ,t)=(\mathcal {L}+\mathcal {L}_\mathrm{qdiff}+\mathcal {L}_\mathrm{mdiff})Q(z,\alpha ,t). \end{equation} \tag{ 13 }$$

The first differential operator on the rhs represents a classical drift

$$\begin{eqnarray} &&\mathcal {L} = \displaystyle {{\rm i}\eta (t)\,{\rm e}^{-{\rm i}\omega _dt}\partial _\alpha +{\rm i}\partial _\alpha \left[\alpha +3g\frac{z+z^*}{1+|z|^2}\right]} \nonumber \\ &&+{\rm i}\partial _z[-\Omega z+g(1-z^2)(\alpha +\alpha ^*)]+{\rm c.c.}, \end{eqnarray} \tag{ 14 }$$

the second quantum diffusion

$$\begin{equation} \mathcal {L}_\mathrm{qdiff}={\rm i}g\partial _\alpha \partial _z(1-z^2)+{\rm c.c.}, \end{equation} \tag{ 15 }$$

and the final term describes the measurement induced diffusion

$$\begin{equation} \mathcal {L}_\mathrm{mdiff}=-\frac{\kappa }{2}(\partial _\alpha -\partial _{\alpha ^*})^2. \end{equation} \tag{ 16 }$$

From Q(z, α, t) we find the field Q-function by integrating out the spin coherent state variable z. The form of the second term $\mathcal {L}_\mathrm{qdiff}$ prevents sub-Planck structures from forming in agreement with the Heisenberg uncertainty relation [51], i.e. the origin of this term is the quantum pressure which becomes very large whenever the distribution is tightly squeezed. This explains the quantum nature of this effect (and this term indeed vanishes in the classical limit ℏ → 0). Note further that $\mathcal {L}_\mathrm{qdiff}$ mixes spin and field variables. The last term $\mathcal {L}_\mathrm{mdiff}$ stems from the noise resulting from the position measurement and only contains field variables. In this model, where $\hat{x}$ is being measured, it has a diffusive effect on its conjugate variable $\hat{p}$. Its form is the same in the Fokker–Planck equation obeyed by the Wigner distribution which explains the smoothening effect seen in figure 2 as κ is increased. To summarize: $\mathcal {L}$ reproduces the classical dynamics, $\mathcal {L}_\mathrm{qdiff}$ gives quantum corrections, e.g. prevents sub-Planck structures from forming and generates qubit–field entanglement, and finally $\mathcal {L}_\mathrm{mdiff}$ describes the back-action of the measurement.

As already mentioned, chaotic evolution typically implies a rapid spreading of the phase space distribution. The 'size' of the distribution can be estimated from the width

$$\begin{equation} \Delta _{xp}=\Delta x\Delta p=\sqrt{(\langle \hat{x}^2\rangle -\langle \hat{x}\rangle ^2)(\langle \hat{p}^2\rangle -\langle \hat{p}\rangle ^2)}, \end{equation} \tag{ 17 }$$

i.e. the product of the quadrature uncertainties. Figures 2 and 3 give the distributions as evolved according to the master equation (6), but as discussed above we may also consider the distributions $Q_i(\tilde{x},\tilde{p},t)$ and W_i(x, p, t) for single trajectories i. Each one of these, corresponding to a single stochastic simulation of the unravelled master equation, gives a width $\delta _{xp}^{(i)}=\sqrt{(\langle \hat{x}^2\rangle _i-\langle \hat{x}\rangle _i^2)(\langle \hat{p}^2\rangle _i-\langle \hat{p}\rangle _i^2)}$ which is used to define the average width

$$\begin{equation} \delta _{xp}=\frac{1}{N}\sum _{i=1}^N\delta _{xp}^{(i)}. \end{equation} \tag{ 18 }$$

Δ_xp and δ_xp represent, respectively, the non-selective width in which measurement results are not recorded and selective width in which the measurements are projective and information about the results of a measurement are fed back into the system. For the same decay rates κ as in figures 2 and 3, we display the time evolution of Δ_xp and δ_xp in figures 5(a) and (b). For a comparison to the semi-classical results we have also calculated Δ_xp within the TWA framework using (8). In addition, in the appendix we have provided snap shots of the individual x- and p-quadrature probability distributions in order to show how they look. Figure 5 reveals numerous interesting phenomena.

• 1.  
  Even for large decoherence rates κ, the position measurement does not localize the state $\hat{\rho }_f(t)$ such that it can be said to approximate the semi-classical result. This is in stark contrast with results found for the Duffing oscillator [33]. In [18] a similar breakdown of the quantum-classical correspondence was conjectured to be due to the quantumness of one of the subsystems (two-level system). Below we will motivate this by calculating the qubit–field entanglement, and it will be shown that despite strong decoherence substantial non-classical correlations between the qubit and the field survive.
• 2.  
  A much more distinct localization effect is found in the selective width δ_xp. For large times, the quantum distributions can even remain more localized than the semi-classical one given by the TWA. We have verified that the localization occurs in both conjugate variables x and p, and not only in the measured position x (see also the appendix).The localization here occurs due to quantum jumps. Every time a jump occurs, the wavefunction ψ(x, t) is multiplied by x. Such a 'measurement' interrupts the evolution and these stochastic kicks of the wavefunction cause a destructive interference leading to the localization effect. This is reminiscent of dynamical localization, which is the analogue in time of 1D Anderson localization in space, and occurs because of random (but coherent) interference that prevents the wavefunction from spreading (the latter has in fact been seen in a periodically 'kicked' optical cavity filled with a nonlinear medium [52]). Here the kicks result from measurements. However, it is interesting to note that the wavefunction of a single simulated trajectory ψ_i(x, t) is always a pure state and so the mechanism for dynamical localization is perhaps quite close even if the origins are different. Meanwhile, the full quantum probability distribution, as given by the solution of the master equation, and which also coincides with the average over all trajectories, does not show the same localization effect. This is not surprising since different trajectories are added incoherently and so one should not expect destructive interference.
• 3.  
  In figure 5(a) we see that at very short times (t ≲ 4) the semi-classical (TWA) and quantum evolution both agree. The quantum and semi-classical curves then diverge, but the various quantum curves continue to agree with each other until the time τ_D ≈ 15 irrespective of κ. This dwell time is a characteristic of open chaotic systems [53, 54] and represents the period during which the open (chaotic) system is transparent to decoherence and evolves similarly to the corresponding closed system. Normally, during the dwell time the system evolves essentially classically and is therefore independent of ℏ [54]. In our case we again find a breakdown from the typical behaviour as the semi-classical broadening diverges from the quantum one long before the dwell time is over.
• 4.  
  At times greater than the dwell time, both figures 5(a) and (b) show that the isolated quantum system (κ = 0) diverges from the open ones (κ ≠ 0). However, the latter show only a very weak dependence on κ (provided that it is not vanishingly small). This is particularly true in the case of figure 5(a). A phenomenon that is perhaps related to this effect is known from studies of Loschmidt echoes in time-reversed open quantum systems where it was discovered experimentally [55], and subsequently understood theoretically [56], that there is a universal regime in which the effects of decoherence are dependent only upon the dynamics of the underlying closed classical system and are, paradoxically, independent of the system-environment coupling strength κ. This regime is known as the Lyapunov regime, the name indicating that decoherence is dependent only upon the Lyapunov exponents characterizing entropy growth in the underlying classical chaotic system [57].

One explanation for the failure of the semi-classical description to reproduce the fully quantum ones even in the presence of decoherence comes from the fact that some quantum properties are maintained despite the decoherence. Since the state $\hat{\rho }(t)$ is in general mixed, standard entanglement measures like the von Neumann entropy cannot be applied to estimate qubit–field entanglement. Neither can the concurrence be utilized since the two subsystems do not constitute two qubits. As a measure of entanglement we use instead the negativity [58]

$$\begin{equation} E_N(\hat{\rho })=-\sum _i\frac{|\lambda _i|-\lambda _i}{2}, \end{equation} \tag{ 19 }$$

where λ_i is the ith eigenvalue of the partial transpose $\hat{\rho }^{{\rm T}_A}$ of $\hat{\rho }$ [59] (here 'A' represents one of the two subsystems—qubit or boson mode). The definition (19) means that the negativity is the sum of negative eigenvalues λ_i. In principle, $E_N(\hat{\rho })<0$ is only a necessary (but not sufficient) condition for entanglement, but, nevertheless, we will take it as an indicator of quantum correlations. We have a lower bound $E_N(\hat{\rho })\ge 1/2$. The negativity for a relatively large κ and for early times is displayed in figure 6(a). At later times, our numerics indicate that $E_N(\hat{\rho })$ shows rapid variations around the value ∼−0.05 and never decays to zero. Thus, the entanglement is not completely destroyed by the decoherence. It has been shown that for closed quantum systems, non-zero entanglement can be a signature of quantum chaos [6, 60]. In general, chaotic dynamics shows large entanglement in closed systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Typically, given a Lindblad term $\mathcal {L}[\hat{A}]$ with $[\hat{A},\hat{H}]=0$ the state $\hat{\rho }(t)$ tends to relax to a steady state which is the ground state of the parent Hamiltonian $\hat{A}^\dagger \hat{A}$ [61]. Here, since $[\hat{x},\hat{H}]\ne 0$ there is an interplay between the two terms, unitary Hamiltonian time evolution and non-unitary decoherent time evolution, and the state is not expected to approach some pure steady state. Indeed, the state $\hat{\rho }(t)$ is highly mixed as can be verified by calculating the purity $P(t)=\mathrm{Tr}\, [\hat{\rho }^2(t)]$ [62]. Although not shown here, we have calculated P(t) for different κs and only in the special case where κ = 0 is the state pure, i.e. P(t) = 1 for κ = 0 and P(t) < 1 otherwise.

As mentioned above, another signature of non-classical features is the negative ratio N₋(t) of the Wigner function [63]. The rule-of-thumb is that decoherence tends to suppress non-classical effects, such as; squeezing, sub-Poissonian photon statistics, and coherent superpositions. We also saw an example of this in figure 2. Examples of the negative fraction N₋(t) are presented in figure 6(b). While N₋(t) is suppressed, there is always a fraction that survives. Compared to the negativity, N₋(t) seems somewhat more sensitive to decoherence.

We have made much in this paper of the discrepancy between the quantum and classical behaviour of our system. However, it is worth noting the following empirical fact: when the dynamics in phase space is averaged over time we find that classical and quantum come back into correspondence. This is illustrated in figure 7 where it can be seen that the total region of phase space explored by the TWA during the dynamics is very similar to the total region of phase space explored by the Wigner distribution, even though, as we saw above in figure 5, at any instant in time their distributions in phase space are quite different. Furthermore, similar looking structures seem to be present in both the quantum and classical distributions. It is noteworthy that for the quantum (Wigner) case in figure 7 we have used unitary time evolution (κ = 0), and so from our previous discussion one might otherwise expect the two distributions to be very different. Instead, the quantum evolution obeys the dynamical limits set in phase space by the classical evolution. In a closed system we might expect this correspondence on the grounds of energy conservation, but here we have an open system which is externally pumped and such arguments do not apply. This correspondence between the time-averaged quantum and classical dynamics has been noted before, e.g. in experiment [6] that reconstructed the time evolution of the Wigner function from measurements. They found that the time-averaged quantum distribution spread throughout phase space but respected the boundaries between the regular and irregular regions with high fidelity.

Zoom In Zoom Out Reset image size

4.2. Detection

We have seen examples of how the effects of decoherence come into play for chaotic quantum dynamics. Therefore, a most relevant question concerns experimental observation of both quantum chaotic dynamics itself, and the effect decoherence has upon it. Thus far we have considered the phase space distributions and especially their widths. The full phase space distribution is in principle attainable experimentally via state tomography [64, 65] (see also [6]). Such homodyne detection relies on experimental repeatability, i.e. repeated measurements on copies of equally prepared systems. Since the system is manifestly chaotic and nonlinear, small fluctuations in the preparation will cause large variations in the evolved state and the measured results may seem random. The measurement result will therefore more likely represent a state averaged over many 'realizations'. This, however, is in itself a signature of the chaotic nature. Thus, state tomography should render a phase space distribution with seemingly random structures.

We may ask whether there are other quantities that are more easily obtained and still show chaotic characteristics? It turns out that even the statistical properties of the intensity of the output cavity field contain information of the quantum chaotic behaviour. The phase space distribution, apart from spreading, will bunch back and forth along the x-direction. When passing the centre at x = 0, population transfer between the two qubit components takes place, i.e. adiabaticity breaks down. For short times, when the wave packet is still localized, this will be manifested as sudden jumps in $\langle \hat{\sigma }_z\rangle$. After the delocalization, the clear transitions between the two components turn into rapid fluctuations. This will also be reflected in large fluctuations in the field intensity $\langle \hat{a}^\dagger \hat{a}\rangle$. Classically, the ergodic evolution implies that a phase space trajectory frequently intersects all accessible 'rings' with radii $r=\sqrt{p^2+x^2}/2$ in phase space. Thus, a classical solution spends some time at all accessible distances from the origin of phase space, and such chaotic motion results in pronounced field intensity fluctuations. For a coherently driven system, the intensity of the output field at time t, $\langle \hat{a}_\mathrm{out}^\dagger \hat{a}_\mathrm{out}\rangle _t$, leaking through the cavity mirror is proportional to the cavity mode intensity $\langle \hat{a}^\dagger \hat{a}\rangle _t$ [66]. Consequently, direct photon detection of the output field gives a handle on the cavity photon fluctuations.

In figure 8 we compare the photon variance of the full field state $\hat{\rho }_f(t)$ (i.e. not for single trajectories) for a constantly pumped system (a), with the one driven with an oscillating amplitude (b). In (a) we show the variance $\Delta _n=\sqrt{\langle (\hat{a}^\dagger \hat{a})^2\rangle -\langle \hat{a}^\dagger \hat{a}\rangle ^2}$ and in (b) the scaled variance $\Delta _n/\langle \hat{a}^\dagger \hat{a}\rangle$. The reason why we do not divide the variance in the regular case (a) is because $\langle \hat{a}^\dagger \hat{a}\rangle$ periodically becomes approximately zero and would cause divergences. For the constant driving case we find a nicely oscillating behaviour, while in the chaotic case the variance shows much greater fluctuations. Another difference is that decoherence enhances field intensity fluctuations in the regular regime, while it suppresses it in the chaotic regime. As a comparison we also display the TWA result. As is seen, without decoherence the quantum variance exceeds the semi-classical one. For strong enough decoherence, the quantum variance may actually be smaller than the semi-classical variance.

Zoom In Zoom Out Reset image size

In this paper we have proposed using a circuit QED setup for studies of quantum chaotic dynamics in a system with no clear classical limit. Special attention was paid to the influence of decoherence deriving from continuous quadrature measurements of the resonator field. Since the measurement back-action only induces decoherence and no dissipation, the corresponding semi-classical equations of motion are transparent to this effect. We demonstrated that the set of semi-classical equations of motion predict classical chaos. In general, decoherence, as it is known to diminish quantum coherence, may explain how classical chaos can appear from quantum mechanics. This is typically the case for systems with a classical limit. However, we showed that this method of achieving classical–quantum correspondence fails in our bipartite model where one of the two subsystems is manifestly quantum. In particular, within a wide range of decoherence rates κ we found that the semi-classical results could not be retrieved (except by averaging over time). By calculating the negativity we concluded that this discrepancy between quantum and semi-classical chaotic evolution could originate from the presence of truly quantum properties such as entanglement even in the presence of strong decoherence. While the full quantum state $\hat{\rho }(t)$ did not show much localization, the states appearing from 'selective measurements' gave a form of localization distinct in origin from standard dynamical localization but actually rather close in mechanism to the stochastic Schrödinger equation we employed to describe the measurements. We further argued that the chaotic nature of the system should show up in the photon statistics which are readily measured in experiments. We note that the amplitude of the information transfer rate κ should be controllable by adjusting the homodyne detection.

It should be pointed out that our results apply to other chaotic models exposed to decoherence. Maybe the most relevant example is the Dicke model [67] describing a set of N two-level qubits coherently interacting with a single cavity mode. The ultrastrong coupling regime was recently achieved experimentally in an effective Dicke model [68], and by exciting the system sufficiently it can be expected to enter a chaotic regime, where our predictions should be detectable in the output cavity field. Our analysis assumes pure decoherence and no dissipation, but we have also tried other forms of the Lindblad terms including, for example, photon decay, and found similar results. When the system is dissipative, one must keep in mind that photon decay itself may cause a localization of the phase space distribution due to energetic considerations that are different from the localization discussed here.

We thank Alexander Altland, Shohini Ghose and Göran Lindblad for enlightening discussions, and anonymous referees for useful comments and suggestions. JL acknowledges support from the Swedish research council (VR), Kungl Vetenskapsakademien (KVA), and Deutscher Akademischer Austausch Dienst (DAAD). DO acknowledges support from NSERC (Canada).

In section 4 in the main text we presented results concerning the area Δ_xp in phase space occupied by the quantum and classical distributions (and also the area δ_xp associated with a sum over single stochastic trajectories). However, since our measurement scheme only measures the x-quadrature it is interesting to know what the typical behaviour of each individual quadrature (the so-called residuals of the phase space probability distribution) is rather than just the combined uncertainty shown in figure 5. For this reason we display in figures A1–A3 snapshots of the probability distributions for the two field quadratures individually at time t = 200. It is interesting to note that there is more fine structure present in the quantum probability distribution for κ = 0 than for κ = 0.05, which is a result of the suppression of quantum interference. In the semi-classical case, the smoothening effect due to the quantum diffusion $\mathcal {L}_\mathrm{qdiff}$ is absent resulting in much finer structures, see figure A3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size