Quantum chaos in an ultra-strongly coupled bosonic junction

The classical and quantum dynamics of two ultra-strongly coupled and weakly nonlinear resonators cannot be explained using the Discrete Nonlinear Schr\"odinger Equation or the Bose-Hubbard model, respectively. Instead, a model beyond the Rotating Wave Approximation must be studied. In the classical limit this model is not integrable and becomes chaotic for a finite window of parameters. For the quantum dimer we find corresponding regions of stability and chaos. The more striking consequence for both classical and quantum chaos is that the tunneling time between the sites becomes unpredictable. These results, including the transition to chaos, can be tested in experiments with superconducting microwave resonators.

The semi-classical and quantum dynamics of two ultra-strongly coupled nonlinear resonators cannot be explained using the Discrete Nonlinear Schrödinger Equation or the Bose-Hubbard model, respectively. Instead, a model beyond the Rotating Wave Approximation must be studied. In the semi-classical limit this model is not integrable and becomes chaotic for a finite window of parameters. For the quantum dimer we find corresponding regions of stability and chaos. The more striking consequence for both semi-classical and quantum chaos is that the tunneling time between the sites becomes unpredictable. These results, including the transition to chaos, can be tested in experiments with superconducting microwave resonators. The unquestionable relevance of the Discrete Nonlinear Schrödinger equation (DNLS) extends beyond the theoretical characterization of waves in nonlinear media [1,2], to describe also a rich variety of physical phenomena ranging from biological physics [3] to gravitational analogues [4]. It is therefore not surprising that the theoretical predictions of the existence of, e.g. discrete solitons [5], vortices [6] and self-trapping [7] can be experimentally tested in a wide variety of setups, such as arrays of optical waveguides [8], polaritons [9] or Bose-Einstein condensates [10]. In this respect, the DNLS equation bridges nonlinear science and quantum many body physics, since it is the semi-classical or many boson limit of the Bose-Hubbard (BH) model, describing amongst other things tunneling of ultra-cold atoms in optical lattices.
The most basic realization of the DNLS is the dimer, formed by two weakly coupled nonlinear resonators. The two-site DNLS is integrable and exhibits a transition [11] from linear oscillatory dynamics (Rabi regime) to self-trapping (localization). Between these limiting cases lays the Josephson regime, where the quantum equivalent of the dimer, the two-site BH model [12] behaves as a bosonic Josephson junction [13]. Furthermore, the signatures of the classical symmetry breaking bifurcation can also be observed in the quantum limit [14]. In that sense the dynamical behavior in the DNLS resembles the Mott-Superfluid transition in the BH model [15].
In this work we study a more general model of two oscillators with amplitude ψ with coupling no longer captured by perturbation theory: the ultra-strongly coupled (USC) bosonic junctioṅ for k = 0, 1. The USC model ( θ = 1 above) becomes the DNLS [θ = 0 in (1)] with nonlinearity strengthγ in the limit of weak couplings |J/ω| 1, by means of the Rotating Wave Approximation (RWA). However, in the ultra-strong coupling regime, where |J| |ω|/2 (see suppl. mat.), the RWA breaks down and new physics is found. The quantum equivalent of (1) are now the ultra-strongly coupled nonlinear resonators, with Fock operators [â k ,â † j ] = δ jk of both oscillators with frequency ω. This quantum dimer lacks a superselection rule and no longer conserves the number of particles, N = kâ † kâ k , just like the classical model (1) no longer conserves the norm N = k |ψ k | 2 , when θ = 1.
The semi-classical and quantum versions of these ultrastrongly coupled bosonic Junctions are of great relevance in the study of superconducting quantum circuits [16]. Labeled as quantum optics on a chip, quantum circuits have reproduced most interesting features of cavity quantum electrodynamics (QED) using photons and superconducting qubits as simulators of light and atoms. One of their greatest advantage is the possibility of pushing the light-matter interaction strength close to the energy of the bare frequency transitions -the ultra-strong coupling. Quite recently this new regime has been experimentally demonstrated in these and other solid state setups [17]. Along this letter the ultrastrong refers to the coupling between two nonlinear bosonic modes. Couplings between superconducting resonators, as the ones we are going to discuss in this letter, recently has been reported experimentally [18].
What Physics can be expected in these new coupling regimes? Our study of the ultra-strongly coupled semiclassical and quantum bosonic dimers reveals that both systems experience a transition to chaos for negative values ofγ. The chaotic regions are finite and can be characterized both spectrally and in phase space, but the most arXiv:1309.0391v2 [quant-ph] 12 Mar 2014 clear signature is the change in the self-trapping dynamics, found in the unpredictability of the tunneling time. These features can be observed using quantum circuits, either in the few-photon or in the semiclassical regime.
The semi-classical limit.-In the limit of large number of excitations,n k := â † kâ k 1, where quantum fluctuations become small, the semiclassical dynamics of the coupled resonators can be approximated using coherent states. Replacing the mean-field valueŝ a k → â k := ψ k in Eq. (2), we obtain a classical Hamiltonian for the dynamical variables ψ k which, by the Hamilton equationsψ k = −i∂ ψ * kH , evolve according to the ultra-strongly coupled bosonic junction equation (1).
In the DNLS limit (θ = 0) the system is integrable and both the energy and the total number of excitations are conserved. Moreover, the equations are symmetric under the transformationγ → −γ, ψ → ψ * and ψ k → ψ k exp(ikπ) (or J → −J). This symmetry and conservation law disappear when we consider the USC model, θ = 1, which is no longer integrable. This has dramatic consequences for the dynamics. We have calculated the semiclassical evolution of the population imbalance ρ(t) := n 0 (t) − n 1 (t) using Eq. (1), comparing the results obtained from the DNLS with the USC model. Along this work we set J = 1 and start with the initial condition n 0 (0) = N 0 , n 1 (0) = 0. This leaves as only free parameters the relative coupling J/ω and interaction strengths γ =γN 0 /J. The normalized minimal imbalance ρ min = min t [ρ(t)/N (t)] is a witness of self-trapping: a value ρ min −1 indicates an oscillating dynamics where particles eventually tunneled to the opposite site, while the imbalance remains locked around ρ min 1, when self-trapping dominates the dynamics.
In Fig. 1(a) we plot ρ min for the USC model, after a sufficiently long time t ∼ 100/J, starting with N 0 = 1. For weak coupling J/ω → 0, the DNLS dimer is recovered and the self-trapping transition happens at the analytic value |γ RWA c | = 4 [11], denoted by white dashed lines in Figs. 1(a-c). For increasing coupling strength, J/ω, and positive γ, self-trapping is observed at slightly smaller nonlinearities γ. More interesting is the behavior for negative γ. Increasing J/ω from zero, the transition is shifted to values |γ| > |γ RWA c |, reaching a minimum at J/ω ≈ 0.1. Around this value we begin to observe strong irregular oscillations of ρ min .
To better understand the dynamics in this parameter region we compute the normalized spectral density defined in terms of f k (ν), the Fourier transform of ψ k (t). Figs. 1(b) and (c) show g(ν; γ) for the RWA case and the full mode respectively, at fixed J/ω = 0.5, with dashed vertical lines delimiting the analytical prediction |γ RWA c | = 4. In Fig. 1(b) we see uniform and well separated lines, indicating that the dynamics is dominated by only a few frequencies, which we identify with the nonlinearly shifted normal modes or the decoupled localized mode. However, in Fig. 1(c) we observe two broad windows with a huge number of frequencies involved in the dynamics, indicating the presence of chaos for negative values of γ.
To get an insight about why chaos can be only found for negative γ, it is worth looking at the linear modes of the systems and their nonlinear continuations. First, we search for linear solutions (γ = 0) of Eq. (1), and we find the corresponding eigenvalues {ν} and eigenvectors of the system (see suppl. mat.). Two modes are symmetric and the other two are anti-symmetric, each having a positive and a negative eigenfrequency. For finite γ, we are able to continue both types of states into the nonlinear regime, finding that Since ν has to be real, in the USC model (θ = 1) there exists a window in parameter space where one of these modes exists: In Fig. 1(c) we also plot the eigenfrequencies given by Eq. (5)  However, for −8 < γ < −4 the antisymmetric modes are missing, and a second region of chaos appears. In Fig.  1(a) and (c) it is possible to observe for γ < −4, first an excitation of the nonlinear self trapped state, and then a weaker chaotic region.
To obtain a more complete characterization of chaos we have computed Poincaré sections in the ρ-φ reduced phase space of the system, where φ(t) := arg[ψ 0 ]−arg[ψ 1 ] is the phase difference between oscillators. We plot points for which N = N coincides with the norm averaged over the integration interval. We start in Figs. 2(a-b) with the DNLS or RWA (θ = 0), for which the norm N (t) = N is conserved, and we have two types stationary orbits depending on the value of the initial imbalance: periodic orbits and self trapped states. The situation is quite different in the USC model, where N is not conserved. For positive values of γ, the Poincaré sections show deformed tori corresponding to the quasiperiodic motion of the oscillators [ Fig. 2(c)], but for negative γ, these tori coexist with chaotic trajectories [ Fig. 2(d)]. The chaotic nature of these orbits was confirmed by the computation of Lyapunov exponents, which were found to be positive for the chaotic trajectories shown in Fig. 2(d).
Quantum dynamics.-A similar study has been done for the quantum models in Eq. (2). Starting with an initial state that corresponds to an imbalanced Fock state |ψ = |n 0 (0) ⊗ |n 1 (0) , with N 0 = n 0 (0) + n 1 (0), we simulated the evolution of the number of photons, n 0,1 (t), and their imbalance, using exact diagonalizations in a truncated Fock basis.
Our first tool for analyzing the dynamics is again the spectral density g(ν; γ) from Eq. (4), which is computed replacing f k (ν; γ) with the Fourier transform of n k (t) for a given γ. The spectral density is shown in Fig. 3 for positive frequencies [19]. In the quantum case with the RWA (θ = 0) we observe the same symmetry γ → −γ [cf. Fig. 3(a)] as in the classical DNLS dynamics. The absence of this symmetry in the USC model (θ = 1) also manifests in Fig. 3(b) with the appearance of a multitude of frequencies for negative γ, in stark resemblance of our semiclassical signature of quantum chaos.
As in the semiclassical limit [cf. Fig. 1] the selftrapping transition for γ < 0 is shifted due to the excitation of chaotic trajectories. To resolve this transition more clearly, we have computed a dimensionless tunneling time τ := min(Jt : ρ = 0) defined as the time at which the population imbalance first changes its sign, for initial conditions n 0 (0) = N 0 , n 1 (0) = 0. Let us first discuss the situation with a large number of particles, N 0 = 17, where the quantum and semiclassical models are expected to converge. The classical self-trapping regions correspond to the shaded areas in Fig. 3(c), and the corresponding semiclassical (quantum) tunneling times are shown in red (black) lines, either dashed (θ = 0, RWA) or solid (θ = 1, USC model). For positive γ we observe that both in the RWA case and the USC model τ grows steeply when entering the shaded region. The only difference is that the quantum model does not lower the self-trapping transition. For negative γ the RWA curve behaves symmetrically, whereas in the USC model we find an irregular behavior. Furthermore, we also observe a small window of regularity in the quantum curve, which corresponds to the same phenomenon of a small intermediate region of quasi periodicity observed in the semiclassical case.
For smaller number of particles, such as N 0 = 2 (green curves), quantum fluctuations become relevant and we enter the so-called Rabi regime of the dimer [13]. The first consequence is that self trapping is possible only for very large onsite interaction γ, outside the range of values considered in this work. The second consequence, the smoothness of the curves, is a clear signature of the lack of chaos. For N 0 = 2 the system is certainly far away from be semiclassical.
To obtain better insight in the tunneling dynamics, we integrated the numerical diagonalization of (2) for t ∈ [0, 1000] and ∆τ for J/ω = 0.5, n 0 (0) = N 0 = 17, γ = −1. The tunneling times ∆τ i = τ i+1 − τ i were defined as the time differences between to consecutive roots of the population imbalance ρ(Jτ i ) = 0. The probability distribution p(∆τ ) is shown in Fig. 4. For the case of RWA, two main peaks can be observed, which are caused by Rabi-oscillations as well as their collapse and revival. The tunneling time distribution for USC is much broader, showing that the tunneling times become irregular, although there seam to be reminiscences of the two former peaks with a main peak for lower tunneling times. This behavior is replicated in the structure of the Fourier transform (cf. Fig (3)), where the ultrastrong coupling exhibits, that much more frequencies are involved in the dynamics.
Discussion.-The experimental achievement of the ultra-strong coupling regime in light-matter interactions is a technological revolution in the field solid state quan-tum optics [17], quantum information [20,21] and quantum simulation [22,23]. However, these experiments have also fundamental theoretical impact, resuscitating old discussions [24], about the integrability of the Rabi model or the emergence of chaos, which see the light with new methods to address them [25].
The present work studies the presence of chaos in the dynamics of the classical and quantum dimer model. Our approach, which is based on the Heisenberg picture of observables or expectation values, is complementary to alternative studies of quantum chaos, such as the energy level spacing statistics, and provides a very natural connection between the quantum and classical worlds. These techniques can be exported to other problems, such as the Rabi model [24,26], which may be regarded as a dimer model in which the nonlinearity of one of the wells drops to zero and the other one is raised to infinity, creating a qubit. This particular model has been shown to be integrable in a novel sense [25], admitting the full classification of energy levels and eigenstates. Compared to the dimer model, it seems that this integrability is enough to support more stable and regular dynamics, and indeed preliminary studies show that the features of quantum chaos in our model are absent in the full Rabi dynamics.
Our works studies the necessary extension of the DNLS beyond the RWA regime. We have reported on a minimalistic system supporting chaos both in the semiclassical and quantum domains.These versions of the ultrastrong bosonic junction map the semiclassical and quantum limits of two ultra-strongly coupled nonlinear resonators. An implementation of this model consists of two superconducting coplanar waveguide resonators with an embedded Josephson junction or a qubit providing nonlinearity [27]. Resonator couplings beyond RWA have been experimentally reported [18] and strongly imbalanced states of the microwave resonators can be engineered at will, using, for example, ancilla qubits [28].
For these models and possible experimental setups our work predicts two main features. The first one is a modification of the transition to self-trapping due to the counterrotating terms. The second feature is the emergence of chaos in the photon number dynamics for attractive bosonic interactions. The physical signature of this novel ultrastrong regime is the unpredictability of the tunneling time. We want to emphasize, that in our quantum calculations we did not have to resort to semiclassical or mixed classical-quantum approximation [26], but have been able to observe chaotic dynamics using a moderate number of excitations, while smooth and purely quantum tunneling prevails in the deep quantum regime. Both results can be tested in a coupled-resonator setup, monitoring the field that leaks from the cavity or installing additional qubits that dispersively probe the electromagnetic field.
Finally, a natural continuation of this work is the search of similiar features in extended coupled cavity ar-rays [29,30], studying the propagation of excitations [31], or the self-trapping dynamics in a many-body setup with dissipation. We believe that our results are just a few first examples of the rich non-perturbative theoretical landscape opened by the new generation of solid-state experiments.

SUPPLEMENT -CONTINUATION OF LINEAR MODES
Starting from eq.(1), we look for the eigenvaluesψ k = λψ k . We split ψ k = a k + ib k and obtain the set or, ordering into real and imaginary part, First, for the linear caseγ = 0, we find the four eigenvalues belonging to a symmetric mode with ψ 1 = ψ 2 and an antisymmetric mode with ψ 1 = −ψ 2 . We now search for their nonlinear continuations. Making the Ansatz |ψ 1 | 2 = |ψ 2 | 2 = N/2, we obtain eqs. (5). At this point we should mention, that the solution of ODE (1) separates no longer into independent stationary mode of the type ψ k (t) = exp(iνt)ψ k with time-independent amplitudes ψ k , but has to be constructed like any solution of an homogeneous ODE.

DIFFERENT PARAMETER REGIMES
In Fig. 5 we show that quantum chaos is not only restricted to the ultra-strong coupling value of J/ω = 0.5, but also can be found at J/ω = 0.1. In such a case, the absolute value of the nonlinearity has to be bigger in order to fulfill classical Eq. (6). Therefore, we use γ = 20 for the simulations shown in Fig. 5. We plot full trajectories of ρ(t) with fixed N 0 = 17 and varying initial imbalance. This is done for positive and negative γ in Fig. 5(a) and (b), respectively. For positive γ, the most localized initial conditions, n 1 (0) = 0, 1 and 3, preserve localization, but become more and more affected by quasi-periodicity. When the initial state is less localized, n 1 (0) = 5, localization breaks down and we observe incoherent oscillations with ρ(t) = 0. For negative γ, Fig. 5(b), the picture is different. The most localized initial condition n 1 (0) = 0 has a tunneling probability higher than the less localized initial state n 1 (0) = 1, which decays much slower to ρ = 0. For further increasing n 1 (0), the tunneling time decreases, recovering for n 1 (0) = 3 a similar value like for n 1 (0) = 0 and than dropping further.