Recurrent delocalization and quasi-equilibration of photons in coupled circuit QED systems

We explore the photon population dynamics in two coupled circuit QED systems. For a sufficiently weak inter-cavity photon hopping, as the photon-cavity coupling increases, the dynamics undergoes double transitions first from a delocalized to a localized phase and then from the localized to another delocalized phase. The latter delocalized phase is distinguished from the former one; instead of oscillating between the two cavities, the photons rapidly quasi-equilibrate over the two cavities. These intrigues are attributed to an interplay between two qualitatively distinctive nonlinear behaviors of the circuit QED systems in the utrastrong coupling regime, whose distinction has been widely overlooked.

We explore the photon population dynamics in two coupled circuit QED systems. For a sufficiently weak inter-cavity photon hopping, as the photon-cavity coupling increases, the dynamics undergoes double transitions first from a delocalized to a localized phase and then from the localized to another delocalized phase. The latter delocalized phase is distinguished from the former one; instead of oscillating between the two cavities, the photons rapidly quasi-equilibrate over the two cavities. These intrigues are attributed to an interplay between two qualitatively distinctive nonlinear behaviors of the circuit QED systems in the utrastrong coupling regime, whose distinction has been widely overlooked.
A single quantum emitter strongly coupled to a quantized electromagnetic field can induce a significant interaction among photons [1], which are usually very weakly interacting. The Jaynes-Cummings (JC) model [2], typically realized in a cavity QED system, provides an intuitive understanding. The nonlinearity of the energy spectrum of the JC model causes an additional energy cost to put an extra photon into the cavity [3], which gives rise to the effective interaction energy of photons. When photons are allowed to hop between nearby cavities forming a so-called JC-Hubbard lattice, a new physics of strongly correlated photons emerges [4][5][6][7][8]. In equilibrium, the JC-Hubbard lattice shows critical behaviors that resemble the physics of Bose-Hubbard model [9,10].
An interesting consequence of the competition between the qubit-cavity coupling and the photon hopping is a selftrapping transition [11], as recently observed in an experiment using two capacitively coupled superconducting transmission lines and transmon qubits [12]. The self-trapping transition in tunnel-coupled quantum systems occurs when an on-site interaction energy becomes so dominant that it prevents quantum tunneling through the tunnel barrier [13][14][15]. Likewise, when the nonlinearity induced by the qubit-cavity coupling exceeds the inter-cavity photon hopping, the photon population dynamics undergoes a sharp transition from a delocalized (tunneling) to a localized (self-trapping) regime [11,12].
Meanwhile, the qubit-cavity coupling that is comparable to a qubit transition frequency or a cavity frequency, the socalled ultrastrong coupling (USC), has been recently achieved in experiments [16][17][18][19][20][21]. In the USC regime, the rotating wave approximation leading to the JC model is not applicable, thus the total excitation number is not conserved [22][23][24][25]. The counter-rotating (CR) terms, which are neglected in the rotating wave approximation, play a crucial role in the physics of strongly correlated photons induced by a light-matter interaction. Equilibrium studies on the JC-Hubbard lattice have shown that the USC leads to an Ising-type quantum phase transition and an exotic phase of light [26][27][28][29].
In this Letter, we explore the dynamics of strongly correlated photons in two coupled circuit QED systems in the USC regime. We examine the phase diagram in the parameter space consisting of the qubit-cavity coupling and the inter-cavity photon hopping. We find that as the photon-cavity coupling increases, the dynamics undergoes double transitions first from a delocalized to a localized phase and then from the localized to another delocalized phase. Moreover, the latter phase is characterized by the quasi-equilibration of photon population, despite that the system is finite and closed. We explain the results based on a competition between two qualitatively distinctive nonlinear behaviors of the circuit QED systems in the USC regime: One nonlinear regime, which is commonly associated with the photon-blockade effect and responsible for the first delocalization-localization transition, has been explored in various contexts in previous works [1,3,30]. However, the other nonlinear regime, responsible for the second localization-delocalization transition and the quasi-equilibration dynamics of photon population, has been widely overlooked so far. Interestingly, the same picture also explains the absence of the photon blockade in single photon transfer dynamics studied previously in Ref. [31]. We note that our findings can be observed by combining existing circuit QED technologies used for the JC dimer experiment [12] and for the realization of USC regime [18,32,33]. Model-Our system is described by the Rabi-dimer model,Ĥ consisting of the two Rabi interaction systemŝ coupled to each other via the photon tunneling with amplitude J. The left (L) and right (R) Rabi interaction systems are assumed to be identical with the cavity frequency ω 0 , the qubit transition frequency Ω, and the qubit-cavity coupling strength g. The coupling strength J is assumed to be sufficiently weak (J ω 0 ) as in common experiments [12]. The operatorâ j describes the field mode of the cavity j, and Pauli operatorŝ σ x,y,z j the qubits. Note that the Rabi Hamiltonian (2) contains the CR terms,â iσ − i +â † iσ + i withσ ± i = (σ x i ± iσ y i )/2, in addition to the JC Hamiltonian. We focus our analysis on the resonant case, ω 0 = Ω, where the effective photon-photon interaction is strongest for given coupling strengths. To investigate the photon-localization dynamics, we suppose that the photons are initially localized in one cavity [11][12][13][14][15]. Specifically, we mainly focus on the case where the initial state is of the particular type |Ψ tot (t = 0) = |n i , ↓ L |0, ↓ R with n i > 10, where |n, σ j denotes a product state of n-photon Fock state and σ =↑, ↓ qubit state in the cavity j = L, R. A few remarks are in order: (i) An initial Fock state is just to simplify the discussion. We have investigated the case of initial coherent states [34,35] and the results are essentially the same (see the inset of Fig. 1) within our parameter regime. (ii) A relatively large number of initial photons (n i > 10 in our simulation) is required because otherwise the localization time is known to become too short (the localized phase disappears) [11].
We describe the photon localization-delocalization transition in terms of the unnormalized photon population imbalance parameter, z(t) = N L (t) −N R (t) withN j =â † jâ j , and its time-averaged value z avg = 1 T T 0 z(t) with the "observation" time T . Note that for g/ω 0 1 the total number of photons N tot (t) = N L (t) +N R (t) can be significantly different from the initial number n i of photons, because the CR terms can generate (or destroy) a considerable number of photons from (or to) vacuum. In such a regime the normalized photon population imbalance, , that has been commonly used to distinguish the localized and delocalized regime in the previous studies [11][12][13][14][15], can severely underestimate the imbalance. In our simulation, we set T = 2 × 10 4 /ω 0 , which is sufficiently long in the parameter regime of interest (to be discussed in greater detail below) to distinguish the localized and delocalized phases.
Results- Figure 1 shows the phase diagram of the photon population dynamics in the g-J space determined by z avg . Note that for the Rabi-dimer both g/ω 0 and J/ω 0 become relevant even in the resonant case ω L(R) = Ω L(R) whereas for the JC-dimer at resonance the dynamics is solely governed by the ratio J/g. The phase diagram exhibits a sharp distinction between the localized regime, where the photons are selftrapped, and the delocalized regime, where photons tunnel between two cavities. More importantly, the phase boundary is nonmonotonic: (i) There exists a critical value J c ≈ 0.03ω 0 of J, above which photons are delocalized for all g. Note that J c ≈ 0.03ω 0 is small enough for the tunneling Hamiltonian (2) to be valid for the cavity-cavity coupling. (ii) For J < J c , as g increases, the system undergoes recurrent transitions, first from a delocalized to localized at g c1 and then from the delocalized to another localized phase at g c2 . The critical value of the first transition scales as g c1 ∼ J √ n i , as already shown in Ref. [11]. The second transition, on the other hand, happens at the critical value, g c2 ∼ J c /J (note that Fig. 1 is in a logarithmic scale), which hardly depends on the initial photon number.
A semiclassical approach has been proposed in Ref. [11] that describes the first transition well. However, it breaks down for larger g and completely misses the recurrent delocalization transition. In passing, Ref. [12] ascribes the first delocalization-localization transition to a classical-to-quantum transition in the sense that the collapse-and-revival emerges in the localized regime when either cavity is initially populated with a coherent state, while the delocalized regime is characterized by classical oscillations.
Even more interesting is the dynamical characteristics of the second delocalized phase, clearly distinguished from the first one. To see this, let us turn to the photon number in each cavity N L(R) (t) and the photon population imbalance z(t) in the USC regime, presented in Fig. 2 (a) and (c). The two cavities share almost the same number of photons at any time t after the transient time dynamics; see the inset of Fig. 2  (a). That is, the fluctuations around z avg = 0 are highly suppressed. It is intriguing to find that the system starting from an imbalanced photon population distribution between the cavities rapidly equilibrates to share an equal number of photons, given that our model includes only two lattice sites without any contact to a bath. The quasi-equilibration of photon dynamics is accompanied by the depolarization of the qubits, σ z L(R) (t) ∼ 0; see Fig. 2 (d). Experimentally, the qubit depolarization is a useful signature of the quasi-equilibration as the qubit state is usually easier to probe than the cavity photon number.
The effect is referred to as quasi-equilibration because a closed quantum system requires a recurrence of dynamics at a finite time τ r [36]. However, it is remarkable that according to our numerical simulation τ r T = 10 4 ω −1 0 . For the circuit QED system, the transmission line resonator and the transmon qubit have a few gigahertz frequency, while relaxation rates are typically in a few megahertz. Therefore, τ r is much longer than any time scales relevant to the experiment, practically equivalent to τ r → ∞. We also note an interesting resemblance to the recently predicted quasi-equilibration between two identical finite quantum systems [37], where a common temperature and small fluctuations around time average of any observables are key features. Our finding also provides a concrete example to the recent discussion of the equilibration in a closed quantum system [37,38].
Discussions-We provide qualitative explanations for the dynamical features of the Rabi-dimer model based on the peculiar properties of the dressed states of single Rabi models.
At g = 0, the eigenstates of the Rabi HamiltonianĤ Rabi [see Eq. (2)] is a product state |n |↑ (↓) of a cavity field mode and a qubit. Namely, the field mode is decoupled from the qubit and its dynamical behavior is harmonic (i.e., linear). As g/ω 0 grows from zero until g/ω 0 1, the JC terms start to take effect and play a dominant role over the CR terms. It leads to the coherent superposition of |n |↑ and |n + 1 |↓ and the corresponding √ n-dependent splitting of eigenvalues.
The nonlinearity in this range of coupling g is thus characterized by the √ n-dependence of the energy levels. The socalled JC nonlinearity induces an effective photon-photon interaction and is known to cause the photon-blockade (PB) effect [11]. In the opposite limit (g/ω 0 → ∞ at resonance and g/ √ ω 0 Ω → ∞ in general), the Ω 2σ z -term is negligible, and the eigenstates ofĤ Rabi have the form |n, ±g/ω 0 |± , where |n, α = e αâ † −α * â |n for a complex number α is a displaced Fock state and |± is aσ x -eigenstate. The field mode is thus linear again and the energy spectrum is harmonic [39,40]. One important difference (compared with the g = 0 case) is that the photon number fluctuations in each eigenstate are huge (eventually diverge with g). As g/ω 0 decreases from the infinity to g/ω 0 1, the Ω 2σ z -term tends to induce transitions between |n, ±g/ω 0 |± , whose transition amplitude is determined by n, −g/ω 0 |n, g/ω 0 ∝ e −2g 2 /ω 2 0 L n (4g 2 /ω 2 0 ) where L n is the n-th Laguerre polynomial. The exponential suppression of the transition amplitude between states, |± , due to a state-dependent displacement of an oscillator, |±g/ω 0 , known as the Franck-Condon (FC) effect [41][42][43][44], leads to an exponential suppression of the energy splitting between |n, ±g/ω 0 |± , and it governs the nonlinearity of the field mode in this range of g. It is stressed that in this range the CR terms play a crucial role and enable the vacuum to "erupt" a large number of photons.
The key observation in the above discussions is that one can expect two qualitatively distinctive nonlinear phases of the Rabi model and a transition between them at g ∼ ω 0 (the actual transition point may vary depending on details, such as J). We will refer to them as the photon-blockade (PB) and photon-eruption (PE) phases, respectively. To examine the distinction more closely, we introduce the photon number variance χ = N 2 − N 2 of the eigenstates and the level spacing variance ζ = 1 where E kk ≡ E k − E k are the differences between eigenenergies E k and E k and K are the number of considered eigenstates (ideally K = ∞). Naturally, ζ characterizes the nonlinearity of the system. One can distinguish the two distinct nonlinear phases as a competition between the nonlinearity ζ and the photon number fluctuation χ. For g/ω ≥ 1, the photon number variance of the eigenstates follows that of the coherent Fock state with a coherent amplitude of g/ω 0 ; that is, χ ∝ (g/ω 0 ) 2 [45]. As expected from the above discussions and illustrated in Fig. 3 (b), χ monotonically increases, whereas ζ first increases but then decreases. For relatively small g/ω 0 , the nonlinearity dominates over the photon number fluctuations. It leads to the effective photon-photon interaction-hence the PB phase. For larger g/ω 0 , on the other hand, the photon eruption due to the CR terms give rise to the large photon number fluctuations, which diminish the effective photon-photon interaction-hence the PE phase.
Now it is fairly straightforward to understand the phase diagram in Fig. 1 and the corresponding photon population imbalance dynamics in Fig. 2 of the Rabi dimer model. For example, consider the case of J/ω 0 = 0.01. In the weakcoupling regime (g/ω 0 1), the system is sufficiently linear and photons oscillate back and forth between the two cavities. As g/ω 0 increases, the JC nonlinearity sets in and when g/ω 0 > 0.1 photons are localized in one cavity due to the PB effect. As the coupling increases further so that g/ω 0 > 1, the FC nonlinearity dominates and the system enters the PE phase. In this regime the photons are delocalized again. Unlike the g/ω 0 1 limit, however, the number of photons in each cavity does not oscillate in time but quasi-equilibrates over the two cavities. We remark that the true g/ω 0 1 limit has not been observed in our actual simulation of the Rabi dimer because of the immense computational cost for large g. In such a limit the system becomes completely linear and should exhibit oscillations at a finite frequency (comparable to J).
One remaining question is how the relatively simple system of a Rabi dimer can have a quasi-equilibration state. To address the issue, we note that the initial localized Fock state |Ψ tot (0) involves a wide range of eigenstates |E tot (labeled by an integer index ) as demonstrated in Fig. 3 (b). Combined with the nonlinearity which makes the energy spectrum E tot highly irregular, it leads to the unusually long recurrence time. Indeed, τ r T in both nonlinear phases (the PB and PE phases); see Fig. 2.
The crucial difference between the two nonlinear phases is the photon number fluctuations: In the PB phase the photonhopping J becomes irrelevant while in the PE phase the enhanced photon number fluctuations enable J to equilibrate photons over the two cavities. To see this, we examine N j (t) : Recall that the spectrum involved in the sum is macroscopically wide and highly irregular. Then the off-diagonal terms with the fast oscillating factor e −iE t tend to cancel each other. At long time scales (ω −1 0 t τ r ), one thus expects Since the system is symmetry under L ↔ R, the eigenstates |E tot are either symmetric or antisymmetric. Therefore, one has E tot |N L |E tot = E tot |N R |E tot and hence N L (t) = N R (t) (i.e., quasi-equilibration) in the PE phase. In the PB phase, the eigenstates |E tot are highly localized on either cavity. The overlap with the initial states Ψ tot (0)|E tot is negligible for those states |E tot localized on the right cavity. It implies that N L (t) ∼ n i and N R (t) ∼ 0 in the PB phase. So far we have ignored the A 2 -term due to the electromagnetic vector potential A of the cavity field [46,47],   3. (a) The level-spacing variance ζ (black filled circles, using the left vertical axis) for the 400 lowest levels and the photon number variance χ (solid lines, using the right vertical axis) for the 20 lowest levels of the Rabi model. For g/ω0 > 1, the photon variance increases with (g/ω0) 2 , as indicated by the dashed line. (b) The overlaps between an initial state |Ψ(t = 0) = |ni, ↓ L|0, ↓ R and eigenstates of Rabi dimer Hamiltonian. with r defined by e 4r = 1 + 4D/ω 0 . But the parameters are renormalized as ω 0 →ω 0 = ω 0 e 2r g →g = ge r and J →J = Je 2r . Therefore the A 2 -term tends to decrease the reduced qubit-cavity couplingg/ω 0 by factor e −r keeping the reduced cavity-cavity couplingJ/ω 0 the same. For the true atom-light coupling, the Thomas-Reiche-Kuhn (TRK) sum rule leads to D/ω 0 > (g/ω 0 ) 2 , and the dynamical features we have discussed above are very difficult to observe experimentally. However, in circuit QED systems the underlying physics of qubit-cavity coupling is different, either the TRK sum rule does not apply or the coupling D has an additional suppression factor [46]. Therefore, the A 2 -term does not affect the relevant parameter region.
Remarks-We have mainly focused on the unitary dynamics of the double Rabi systems. In realistic experiments, it is expected that the cavity damping time τ κ (devided by N ) or the qubit decoherence time τ φ sets the observation time T for the recurrent delocalization and quasi-equilibration transition. Note that the recent state-of-the-art experiments [12,[17][18][19]32] have realized τ κ , τ φ > 10 4 /ω 0 . We have briefly examined the cavity damping effect based on the quantum jump approach [48] and taking into account the nontrivial interplay between the damping 1/τ γ and the strong coupling g. As illustrated in the inset of Fig. 1, the localized phase becomes less prominent (as noted in Ref. [11]) but our main result survives small damping. We leave open further extensive studies of the damping effects on the photon localization-delocalization.
We have explored the photon population dynamics in a system of two coupled circuit QED systems in the USC regime. For g ω 0 , the recently observed localized photon dynamics [12] gives way to quasi-equilibration dynamics. It reveals a new qualitatively distinctive nonlinear behavior of the circuit QED system. MJH was supported by the ERC Synergy grant BioQ and the EU STREP DIADEMS and EQUAM and the IT R&D program of MOTIE/KEIT [Grant No. 10043464(2012)]. The Python package QuTiP 2 [49] has been used for simulations. MSK thanks the UK EPSRC. MSC acknowledges the support by the NRF of Korea (Grant No. 2015-003689).