Fock space localization of polaritons in the Jaynes-Cummings dimer model
FFock space localization of polaritons in the Jaynes-Cummings dimer model
Hassan Shapourian and Darius Sadri † Department of Physics, University of Illinois at Urbana-Champaign, Urbana Illinois 61801, USA Department of Electrical Engineering, Princeton University, Princeton New Jersey 08544, USA (Dated: January 25, 2016)We present a new method to study the semiclassical dynamics of the Jaynes-Cummings dimermodel, describing two coupled cavities each containing a two-level system (qubit). We develop aFock space WKB approach in the polariton basis where each site is treated exactly while the intersitepolariton hopping is treated semiclassically. We show that the self-trapped states can be viewed asFock space localized states. We find that this picture yields the correct critical value of interactionstrength at which the delocalization-localization transition occurs. Moreover, the validity of ourWKB approach is supported by showing that the quantum spectrum can be derived from a setof Bohr-Sommerfeld quantization conditions and by confirming that the quantum eigenstates areconsistent with the classical orbital motion in the polariton band picture. The underlying idea ofour method is quite general and can be applied to other interacting spin-boson models.
I. INTRODUCTION
Interacting light-matter systems have opened new pos-sibilities in simulating various types of strongly correlatedphases ranging from the superfluid and Mott insulatingstates to the quantum Hall fluids . The prop-erty that makes the coupled light-matter systems uniquecompared to conventional condensed matter systems isthat they are equipped with optical probes to manipu-late arbitrary states, beyond the ground state, as wellas to monitor non-equilibrium dynamics. The stronglycorrelated light-matter physics have been experimentallyrealized in various implementations of cavity QED withatoms , excitons , or superconducting qubits .The building block of all these systems is the celebratedJaynes-Cummings (JC) model, H JC = ν q σ z + ν c a † a + g ( σ + a + a † σ − ) (1.1)which describes a two-level system (qubit) interactingwith a single mode of the electromagnetic field insidea cavity. The qubit is characterized by the resonancefrequency ν q , and the pseudo-spin operators σ z , and σ ± = ( σ x ± iσ y )/
2. The cavity field with frequency ν c is described by the annihilation and creation operators a and a † which are linearly coupled to the qubit via thecoupling constant g . Throughout this paper, we focuson the resonant case where ν q = ν c = ν . The cavity-qubitcoupling induces an anharmonicity in the spectrum of theJC Hamiltonian that can be qualitatively regarded as aneffective interaction between photons. A tunnel-coupledarray of JC sites would then realize an ideal setup forstrongly correlated photons.The Jaynes-Cummings dimer model (JCDM) is thesimplest, yet non-trivial, system to study interacting pho-tons . The model Hamiltonian H = H JC,R + H JC,L − J ( a † R a L + a † L a R ) , (1.2)describes two identical JC sites coupled through a tun-neling (kinetic) term in the photonic channel a and a † where the subscript L ( R ) specifies the left(right) sites.As the cavity-qubit interaction is increased, the JCDMdisplays a transition from a Josephson oscillating (delo-calized) regime where photons coherently tunnel betweencavities to a self-trapped (localized) regime where pho-tons are frozen inside one cavity. A similar transitioncaused by the Hubbard interaction has been studied ina Bose-Einstein condensate double-well (BECDW) sys-tem . A theoretical study of this transition in JCDMwas originally done in the pioneering work of Schmidt, et al. where a semiclassical picture as well as numer-ical exact solutions were provided. They also proposeda superconducting circuit implementation of the JCDM.This proposal was subsequently realized experimentallyand the transition was successfully observed .In this article, we revisit the localization-delocalizationtransition from a different point of view beyond just nu-merics. We note that the JCDM conserves the total po-lariton number N = n L + n R where n s = σ + s σ − s + a † s a s for s = L / R , left/right JC islands; i.e. [ N, H ] =
0. We usethis property to develop a new approach in the polari-ton basis where each JC site is treated exactly while thehopping of polaritons between sites is treated semiclas-sically. To this end, we introduce a Fock space picturein which the tunneling terms act like hopping operatorsand map the Hamiltonian onto a one-dimensional tight-binding model. The Schr¨odinger equation in this basisobeys a discrete form and hence can be approximatelysolved by a discrete version of the Wentzel-Kramers-Brillouin (WKB) approach . Using our Fock space rep-resentation, we show that the self-trapping transition inthe JCDM can be viewed as a localization transition inFock space . This phenomenon has also been discussedin the context of BECDW . In fact, we find that theWKB approximation maps the JCDM into a four-bandone-dimensional tight-binding Hamiltonian. Each bandlooks like a BECDW with a more complicated interac-tion. These bands are associated with the four possiblepolariton states of two JC sites: upper-upper, upper-lower, lower-upper, and lower-lower polaritons. As wewill see the upper-upper and lower-lower polariton bands a r X i v : . [ qu a n t - ph ] J a n are equivalent to BECDW with attractive and repulsiveinteractions, respectively.Our article is organized as follows. In Sec. II, weintroduce the Fock space picture and find the wave-functions in the WKB limit. In Sec. III, we derive aset of Bohr-Sommerfeld quantization conditions in theform of S ( E n ) = ∮ p ⋅ dq = πh ( n + / ) which yields aquantized energy E n in each polariton band. The mainresults of the paper are in Sec. IV where we study thedelocalization-localization transition using the presentedWKB picture. Finally, we discuss the conclusions andpossible applications of our methods in Sec. V. II. FOCK SPACE REPRESENTATION ANDWKB ANALYSIS
In this section, we write the JCDM in the Fock spaceof photons and show that the Hamiltonian is tridiago-nal. Thus, it can interpreted as a one-dimensional tight-binding model with nonuniform nearest neighbor hop-pings. This representation of the Hamiltonian is going tobe the basis for the exact diagonalization and the WKBanalysis in the upcoming sections.As mentioned earlier, the total polariton number op-erator N = n L + n R commutes with the Hamiltonian andwe can study this model in a subspace with a fixed N .The simplest basis to span this subspace is a set of 4 N or-thonormal states denoted by ∣ n c,L , m L ⟩⊗∣ n c,R , m R ⟩ whichrefers to a state with n c,L ( R ) = a † L ( R ) a L ( R ) photons and m L ( R ) =↑ / ↓ , up/down qubit on the left (right) site. Theonly constraint is that the total polariton number mustbe N. Consequently, any generic state can be decomposedinto a superposition of these states ∣ Ψ ⟩ = ∑ A n c,L m L n c,R m R Ψ ∣ n c,L , m L ⟩ ⊗ ∣ n c,R , m R ⟩ . The polariton imbalance of a state ∣ Ψ ⟩ is defined as theexpectation value of Z = ⟨ Ψ ∣( a † L a L + σ + L σ − L ) − ( a † R a R + σ + R σ − R )∣ Ψ ⟩ (2.1)and the normalized imbalance is defined by x = Z / N .Note that the integer Z ranges from − N to + N in incre-ments of 2. In our WKB analysis the imbalance Z willplay the role of the position. We chose our basis in termsof the eigenstates of imbalance operator ∣ Z, m L , m R ⟩ ,where ∣ Z, ↓ , ↓⟩ = ∣( N + Z )/ , ↓⟩ ⊗ ∣( N − Z )/ , ↓⟩∣ Z, ↓ , ↑⟩ = ∣( N + Z )/ , ↓⟩ ⊗ ∣( N − Z )/ − , ↑⟩∣ Z, ↑ , ↓⟩ = ∣( N + Z )/ − , ↑⟩ ⊗ ∣( N − Z )/ , ↓⟩∣ Z, ↑ , ↑⟩ = ∣( N + Z )/ − , ↑⟩ ⊗ ∣( N − Z )/ − , ↑⟩ (2.2)provided that ∣ Z ∣ < N . Note that for Z = ± N there areonly two states. In the absence of cavity-qubit coupling g , the states ∣ Z, m L , m R ⟩ are fourfold degenerate. This degeneracy is lifted by g through the following term H Z = g ∑ m R √ N + Z (∣ Z, ↑ , m R ⟩⟨ Z, ↓ , m R ∣ + H.c. )+ g ∑ m L √ N − Z (∣ Z, m L , ↑⟩⟨ Z, m L , ↓ ∣ + H.c. ) . (2.3)Therefore, for any fixed N the Hamiltonian in this basiscan be written as H = − J ∑ Z,m L ,m R ( T m L m R Z,Z + ∣ Z, m L , m R ⟩⟨ Z + , m L , m R ∣ + h.c. )+ ∑ Z H Z (2.4) where the hopping terms are T ↓↓ Z,Z + = [( N + Z + )( N − Z )] / ,T ↓↑ Z,Z + = [( N + Z + )( N − Z − )] / ,T ↑↓ Z,Z + = [( N + Z )( N − Z )] / ,T ↑↑ Z,Z + = [( N + Z )( N − Z − )] / , and the first two terms in Eq. (1.1) are dropped as theirsum is constant in the resonant case ν c = ν q . We representwave-functions in terms of a four-component “position”-dependent vector C ( Z ) = [ C ↓↓ , C ↓↑ , C ↑↓ , C ↑↑ ] T , ∣ Ψ ⟩ = ∑ Z,m L ,m R C mLmR ( Z ) ∣ Z, m L , m R ⟩ (2.5)and diagonalize the Hamiltonian in this basis. The spec-tral map is then constructed by tracing out the spin de-grees of freedom as illustrated in Fig. 1; i.e. we compute ∣ Ψ ( Z )∣ = ∑ i ∣ C i ( Z )∣ for every eigenstate and show it asa false color map. The white regions correspond to theclassically forbidden regime where the eigenstates haveexponentially small probability. The localization occurswhen the semiclassical description corresponds to a par-ticle in a double well with the two wells separated bya barrier that the particle has to tunnel through. Thestatement of localization is that the tunneling time fromone half to the other is exponentially large. For instance,for a small ratio of couplings g / J as in Fig. 1(top) thereis no barrier and the localization does not occur, whereasfor larger couplings as in Fig. 1(bottom) a barrier (whiteregions) emerges over some energy ranges, and localiza-tion is possible. The “two wells” correspond to the pho-tons being mostly on one of the two JC sites. When thebarrier is present, the tunneling time to the other siteis exponentially large in N . The spectral map consistsof four bands and each band is associated with a certainpolariton configuration. We shall explain this in more de-tail after deriving the WKB solution. The delocalization-localization transition of the BECDW problem has alsobeen studied in the Fock space representation . Theimportant difference here is that the BECDW Hamilto-nian can be formulated fully in terms of the imbalance FIG. 1. (Color online) Eigenstates and spectrum on the JCdimer model. Color represents the amplitude squared of theeigenstates in Z -space defined by Eq. (2.5). The horizontalaxis is normalized imbalance x = Z / N , while the vertical axisis a scaled eigenenergy. Here g / J √ N = . . N = and the relative phase between two condensates whereasour JCDM has two qubits as extra degrees of freedomthat ultimately lead to four bands in the Fock space rep-resentation.An important remark is that the above Hamiltonian,Eq. (2.4), possesses a Z (left-right) parity symmetrywhose operator P is defined as P∣ Z ⟩ = ∣ − Z ⟩ ; thus, theeigenstates must respect this symmetry, too. This meansthat strictly speaking there is no localized eigenstate andthe localization is purely a dynamical effect in the fol-lowing sense: If the system is prepared with some initialnon-zero imbalance the dynamics would still preserve thisnon-zero imbalance for a long time. The time scale atwhich the imbalance changes sign is divergent as a func-tion of system size.Let us now work out a semiclassical framework to studythe Hamiltonian in Eq. (2.4) in the large-N limit. To this end, we divide the Schr¨odinger equation by NiN ∂ ∣ Ψ ⟩ ∂t = H ∣ Ψ ⟩ , and bring the 1 / N factor into the definition of H . Usingthe basis introduced in Eq. (2.5), we define the effectivePlanck’s constant h = / N and the position (normalizedimbalance) x = Z / N and rewrite the Schr¨odinger equa-tion as ih ∂∂t Ψ ( x, t ) = W ( x ) Ψ ( x ) − B ( x + h ) Ψ ( x + h )− B ( x − h ) Ψ ( x − h ) . (2.6)where Ψ ( x, t ) is the continuum version of the four compo-nent vector C ( Z ) (Eq. (2.5)) and we introduce matrices B ( x ) = J ([( + x + h )( − x + h )] / , [( + x + h )( − x − h )] / , [( + x − h )( − x + h )] / , [( + x − h )( − x − h )] / ) ,W ( x ) = g ′ √ + x σ x ⊗ I + g ′ √ − x I ⊗ σ x . (2.7)where I is a 2 × g ′ = g √ N . Let us introduce the conjugate op-erator to position ˆ p such that [ ˆ p, x ] = − ih . So, theSchr¨odinger equation can be recast into its continuumversion as ih ∂∂t Ψ ( x, t ) = W ( x ) Ψ ( x, t ) − e i ˆ p B ( x ) e i ˆ p Ψ ( x, t )− e − i ˆ p B ( x ) e − i ˆ p Ψ ( x, t ) . The next step is to find a semiclassical ansatz for thewave-function. A comprehensive discussion of solving thediscrete Schr¨odinger equation in the semiclassical WKBlimit can be found in . A similar many-body WKB ap-proach has been applied to the BECDW problem . Abrief review of this procedure is provided in Appendix C.Before we apply the WKB approach to the JCDM, itis essential to write down the matrices in the polaritonbasis given by the modes ∣ Ψ ⟩ = ∣ N ( + x )/ , +⟩ ⊗ ∣ N ( − x )/ , +⟩ , ∣ Ψ ⟩ = ∣ N ( + x )/ , +⟩ ⊗ ∣ N ( − x )/ , −⟩ , ∣ Ψ ⟩ = ∣ N ( + x )/ , −⟩ ⊗ ∣ N ( − x )/ , +⟩ , ∣ Ψ ⟩ = ∣ N ( + x )/ , −⟩ ⊗ ∣ N ( − x )/ , −⟩ , where ∣ n, ±⟩ = (∣ n, ↓⟩ ± ∣ n − , ↑⟩)/√ H JC ∣ n, ±⟩ = ( ν n ± g √ n )∣ n, ±⟩ . (2.8)In the rest of this paper, we refer to the states ∣ Ψ ⟩ , . . . , ∣ Ψ ⟩ by calling them upper-upper (first), upper-lower (second), lower-upper (third), and lower-lower (a) ǫ x − − (2 g ′ + J ) − g ′ √ − JJg ′ √ g ′ + J ) (b) ǫ x − − (2 g ′ + J ) − g ′ √ − JJg ′ √ g ′ + J ) FIG. 2. (Color online) Classically allowed regions (whitearea) of different polariton components for J / g ′ = / (fourth) polariton states or by noting their respective po-lariton index shown inside parentheses. Using this basis,the JC interaction W ( x ) is diagonal W ( x ) = g ′ diag [ A p + A m , A p − A m , − A p + A m , − A p − A m ] where A m = √ − x and A p = √ + x . A perturbativeexpansion of B ( x ) in powers of h is B ( x ) = B ( ) ( x ) + hB ( ) ( x ) + . . . where B ( ) = J √ − x I (2.9) and B ( ) ( x ) = J √ − x ⎛⎜⎜⎜⎝ x +
12 1 − x x + − x − x x + − x x + ⎞⎟⎟⎟⎠ . (2.10)Note that in contrast to the BECDW problem, we havefour polariton modes. The wave-function takes the formΨ ( x, t ) = ⎛⎜⎜⎜⎜⎝ α ( x, t ) e ih S ( x,t ) α ( x, t ) e ih S ( x,t ) α ( x, t ) e ih S ( x,t ) α ( x, t ) e ih S ( x,t ) ⎞⎟⎟⎟⎟⎠ for the Schr¨odinger equation H Ψ = ih ∂∂t Ψ. The real andimaginary parts of the Schr¨odinger equation are straight-forward to compute. Since W ( x ) and B ( ) ( x ) are diag-onal, we get four decoupled equations at zeroth order in h , − ∂ t S i = − B ( ) ii cos ( ∂ x S i ) + W ii ( x ) . From the zeroth-order equations, we can define classicalHamiltonians associated with each polariton mode H = g ′ (√ − x + √ + x ) − J √ − x cos 2 ϕ ,H = g ′ (−√ − x + √ + x ) − J √ − x cos 2 ϕ ,H = g ′ (√ − x − √ + x ) − J √ − x cos 2 ϕ ,H = g ′ (−√ − x − √ + x ) − J √ − x cos 2 ϕ , where at the classical level in each mode ϕ i = ∂S ( ) i ∂x and x are canonical variables such that their Poisson bracketis { x, ϕ } =
1. The four phase variables introduced aboveemerge from the original angle variables associated withqubit degrees of freedom; however, the angles definingqubit states (say in the Bloch representation) are cou-pled to each other and have enormous quantum fluctua-tions, while the emergent angle variables in the presentpolariton basis are decoupled from each other in the clas-sical limit and have exponentially small (in system size)fluctuations in the strong coupling (large g ) limit. Weintroduce the classical velocity, v i = ∂H i ∂ϕ i = B ( ) ii ( x ) sin 2 ϕ i . (2.11)For each Hamiltonian, one can define the allowed regionfor classical solution bound by two turning points wherethe velocity vanishes at 2 ϕ i = π with V li and V hi : V l ( h ) ( x ) = g ′ (√ − x + √ + x ) ± J √ − x V l ( h ) ( x ) = g ′ (−√ − x + √ + x ) ± J √ − x V l ( h ) ( x ) = g ′ (√ − x − √ + x ) ± J √ − x V l ( h ) ( x ) = g ′ (−√ − x − √ + x ) ± J √ − x (2.12)and the classical motion is constrained to be within thehard barriers defined by the above expressions. Figure 2illustrates the allowed region for each polariton mode. Itis helpful to view the JCDM in this basis as a quantumparticle confined inside the classically allowed region. Infact, the bounded regions in Fig. 1 coincide with theclassical potentials in Fig. 2 (see also an illustration ofwave-functions in Fig. 3). Note that in the weakly inter-acting limit (Fig. 2(a)), all bands overlap and there is nogap (forbidden region) in the middle (i.e. no localization)while a forbidden region appears in the strong coupling(Fig. 2(b)).We look for the stationary solution in form of ψ i = α i ( x ) e − ih (cid:15)t e ih S i ( x ) where ψ i is the i -th polariton compo-nent of Ψ ( x, t ) . Note that (cid:15) = E / N is the normalizedenergy where the original eigenvalues of the Hamiltonianis denoted by E . In the parameter regime J / g ′ ≤ − √ (cid:15) can be found easily (see Fig. 2(b)):(i) g ′ √ < (cid:15) < ( g ′ + J ) Only the first polariton component is allowed.(ii) − g ′ √ < (cid:15) < − J and J < (cid:15) < g ′ √ x . In other words, they are notsimultaneously non-zero over the same range of x .(iii) ∣ (cid:15) ∣ < J Both second and third polariton components areallowed. There exists a common range of x in whichboth components are non-zero.(iv) −( g ′ + J ) < (cid:15) < − g ′ √ i − th polariton mode in theclassically allowed regions is ϕ i = ∂ x S ( ) i =
12 Arccos ⎛⎝ (cid:15) − W ii ( x )− B ( ) ii ( x ) ⎞⎠ , (2.13)and in the forbidden region is φ i = ∂ x Q ( ) i =
12 cosh − ⎛⎝ (cid:15) − W ii ( x )− B ( ) ii ( x ) ⎞⎠ , (2.14)where we define S i ( x ) = iQ i ( x ) to get an exponentiallygrowing/decaying wave-function. We note that the solu-tion is oscillating in a classically allowed region whereasit is exponentially decaying in a forbidden region.The first-order correction mixes different modes. For ageneric case, some modes are forbidden and some modesare allowed, so one can write ∂ x S ( ) i = ∑ j B ( ) ij ( x ) α j ( x ) e ih ( S j − S i ) cos ( ∂ x S ( ) j ( x )) B ( ) ii ( x ) α i ( x ) sin ( ∂ x S ( ) i ( x )) (2.15) It is important to keep in mind that for the classicallyforbidden solutions S ( ) is pure imaginary and sinusoidalfunctions must be replaced by the hyperbolic functions. III. BOUNDARY MATCHING ANDQUANTIZATION RULES
In this section, we obtain quantization rules for the en-tire energy spectrum in the strong coupling regime. Thisis to illustrate the usefulness of the semiclassical polari-ton band picture to describe the quantum dynamics interms of classical orbits. For any given energy, one com-ponent is large and others are exponentially small; exceptfor the middle band where both second and third com-ponents are non-zero. As Fig. 4 shows, the quantizationrules agree with the quantum spectrum in all regions aslong as bands do not overlap or J / g ′ < − √
2. A detaileddiscussion of the validity of the WKB approximation isprovided in Appendix G. In short, the classical motionin polariton bands is valid for a wide range of photonnumbers up to such small numbers as N = − J − g ′ ≤ (cid:15) ≤ − g ′ √
2. This band is similar to BECDWwith repulsive interaction. As shown in Fig. 3(a), thereare three regions in this band: (i) (cid:15) < J − g ′ , where the de-localized wave-functions center around x = (cid:15) > J − g ′ , where the localized wave-functions form,in which the probability density is maximum near theboundary points x = ± (cid:15) c, = J − g ′ where the Schr¨odinger equation will not be linear and asa result the quantization rule will be different from theusual Bohr-Sommerfeld formulas. Let us now discuss thequantization condition in these three regions one by one.(i) Delocalized states ( (cid:15) < J − g ′ ): The particle is con-fined in a two-sided potential and the probability densityis large close to the center ( x = ∣ x ∣ < z l where ± z l are turningpoints (cid:15) = V l ( x = ± z l ) , the quantization condition in this (a) ǫ x − − g ′ − J − g ′ √ − g ′ + J −0.200.2 Ψ ( x ) −1 −0.5 0 0.5 1−0.200.2 x Ψ ( x ) (b)(c) FIG. 3. (Color online) (a) Classically allowed region for thelower-lower polariton band. It is divided into the localized anddelocalized regions. The critical energy (cid:15) c, = − g ′ + J (dash-dotted line) separates these two regions. (b) Delocalized state (cid:15) < (cid:15) c, , and (c) localized state (cid:15) > (cid:15) c, . The blue is exactdiagonalization results and green is the WKB wave-function.Here J / g ′ = / N = case is derived to be1 h ∆ S ( (cid:15) ) = ( n + ) π (3.1)where ∆ S = ∫ z l − z l ϕ ( x ) is the classical WKB phase and ϕ ( x ) is the canonical momentum in the fourth bandgiven by Eq. (2.13). As expected, this result is similar tothe familiar Bohr-Sommerfeld quantization for a particleconfined in a potential well . Note that other modescontribute through Eq. (2.14) which are exponentiallysmall and are hence neglected.(ii) Localized states ( J − g ′ < (cid:15) < − g ′ √ x = ±
1. Hence, there are two sets of turningpoints z l and z h : One is bouncing off the ϕ = (cid:15) = V l ( x = z l ) , and −2 −1 0 1 200.01 E rr o r i n W K B ǫ /g ′ −1.76 −1.7400.050.1 E rr o r i n W K B ǫ /g ′ (b)(a) FIG. 4. (Color online) (a) Error in the Bohr-Sommerfeldquantization rules over the entire quantum spectrum. Dif-ferent symbols represent different bands: blue circles areEq. (3.1) and its analog for the upper-upper band, green cir-cles is Eq. (3.2) and its analog for the upper-upper band,and red and black triangles are middle bands. The verti-cal lines represent the following energy scales: red dashed, (cid:15) = ±( g ′ − J ) , black dotted, (cid:15) = ± g ′ √
2, and green dash-dotted, (cid:15) = ± J . (b) Zoomed figure around critical energies (cid:15) c, = − g ′ + J (left) and (cid:15) c, = g ′ − J (right). The red starsare the modified quantization rules at the critical level givenby Eq. (3.4) and its analog for the upper-upper band, theother symbols are the same as in (a). Here J / g ′ = / N = the other is bouncing off the 2 ϕ = π potential, whichis the solution of (cid:15) = V h ( x = z h ) . Provided that theclassically allowed region is z h < ∣ x ∣ < z l , the quantizationrule is found to be1 h ∆ S ( (cid:15) ) = nπ ± e − h ∆ Q ( (cid:15) ) . (3.2)where ∆ S is to be integrated over the classically allowedregion and ∆ Q is the WKB tunneling amplitude (seeAppendix E for explicit expressions). Note that this con-dition is similar to the usual Bohr-Sommerfeld except forthe tunneling term. The important message here is thatthe tunneling rate is related to the energy splitting in theenergy spectrum through δ(cid:15) (split) = e − h ∆ Q ( (cid:15) ) ∣ d ∆ S d(cid:15) ∣ − . (3.3)which is consistent with the quantum spectrum; i.e. forenergies in the range ∣ (cid:15) ∣ < g ′ − J , there are pairs ofeigenvalues that come exponentially close together as thesystem size N = / h is increased. In the classical limit N → ∞ , tunneling is suppressed and we get pairs of fullylocalized degenerate states. Figure 3 illustrates a fewtypical examples that WKB solutions match well withthe exact wave-functions.(iii) Critical region ( (cid:15) ∼ J − g ′ ): Consider a small de-viation λ ≪ (cid:15) = J − g ′ + λ / N .The quantization formula is derived to bearg [ e − ih ∆ S √ π Γ ( / − iχ ) − e πχ / ] = nπ + π / S is calculated over the classical region, Γ ( x ) denotes the Gamma function and χ = λ /√ J ( g ′ − J ) .Notice that the quantization condition in this case ismore complicated than the standard Bohr-Sommerfeldquantization (Eq. (3.1)) and as we show later, this leadsto smaller level spacing and a larger density of statesnear the critical levels compared to regular regions inthe spectrum. This new form of the quantization ruleclose to energy levels dividing localized and delocalizedstates was originally found for the BECDW problem inRef. . As shown therein, an immediate implication ofthis formula is that the quantum break times, at whichthe classical and quantum solutions start to differ, nearcritical regions grow logarithmically with N instead ofalgebraically. This conclusion also applies to our case attwo critical regions of the spectrum for upper-upper andlower-lower bands.Figure 4 summarizes the main results of this section.We plug in the quantum eigenvalues of the Hamltonian tothe quantization condition for the appropriate polaritonband determined by the location of the eigenvalue andsubtract the integer multiple of π from it. The differenceis shown as the defect of WKB quantization. It is evidentin Fig. 4(a) that the WKB error is quite small except forthe critical energy levels (cid:15) c = ±( g ′ − J ) . Close to crit-ical levels, we must use the modified Bohr-Sommerfeldquantization rules, Eq. (3.4), as shown in Fig. 4(b). InAppendix G, we investigate the validity of the polaritonband picture by introducing a quantum-to-classical cor-respondence based on the Husimi distribution. Using thephase space distribution to extract the classical aspectsof the Hamiltonian eigenstates, we show that the proper-ties associated with polariton bands remain valid up tosuch small total polariton numbers as N =
6. We alsonote that the quantization conditions start to fail as wego to a weak coupling regime where the polariton bandsstart to overlap and higher order corrections become im-portant. Nonetheless, the quantization conditions stayvalid away from the overlapping regions in the spectrum(see Fig. S6 in Appendix G).
IV. LOCALIZATION TRANSITION
Before we discuss the localization transition in the po-lariton band picture, let us make few remarks about thetransition in the classical limit of the JCDM. In Ap-pendix A, using a coherent-state path integral formalism V l ( x ) x − g ′ √ FIG. 5. (Color online) Evolution of the minimum potentialcurve V l ( x ) of upper-upper band, Eq. (2.12), as J / g ′ is tunedup from top to bottom as shown by the arrow. The maximumpotential curve V h ( x ) of the lower-lower band is just the mir-ror image of these curves with respect to the horizontal line (cid:15) =
0. The red (dash-dotted) curve for J / g ′ = − √ x = x = ± in the large photon number and large spin limit, we showthat the classical dynamics of the JCDM is described inan eight dimensional phase space, two for photons andtwo for spins per JC site. The resulting equations ofmotion are found to be the same as the factorizationof Heisenberg equations. The critical value of coupling g c , above which localization occurs, is found to be de-pendent on initial spin configurations and its minimumvalue is analytically derived to be g c = J √ N . This isthe same as the quantum mechanical value for g c foundnumerically and analytically (as shown below). Thisresult is in contrast with the previous numerical analysisof , g c ≈ . J √ N , where dynamics in a reduced fourdimensional phase space was only studied. In the samereference it was argued that the difference between classi-cal and quantum values of g c is related to large quantumfluctuations. However, we see that the difference is an ar-tifact of only considering part of the classical phase space,and taking into account the full phase space resolves thisissue.Now, let us turn to the polariton band picture. Aswe have seen in the previous section, the localized statesappear at energies close to where the upper-upper or thelower-lower polariton bands meet the middle bands. AsFig. (5) suggests, the localization can be understood bystudying the changes in the curvature of the V l ( x ) or V h ( x ) functions. Here, we investigate the localizationtransition for the upper-upper polariton band in this sec-tion. The result is exactly the same for the lower-lowerband. For the purely classical solution, where the parti-cle is bound to be in the allowed regions, one can easilyobtain the period of motion for a given trajectory T =∮ dx / v and the average imbalance ⟨ x ⟩ = ( / T ) ∮ xdx / v where the velocity v is given by Eq. (2.11) in terms of theinitial position x where (cid:15) = V l ( x ) . Figure 6(top) shows FIG. 6. (Color online) Top: the classical average imbalance(color code) as a function of J / g ′ for various initial conditions( (cid:15) = V l ( x ) ). Bottom: the expectation value of the imbalance(Eq. (2.1)) as the color code for the eigenstates of the Hamil-tonian close to the touching point between upper-upper andmiddle bands. The dashed line is the classical phase boundarygiven by Eq. (4.1). the average imbalance for various values of the J / g ′ . Thephase boundary is given by gJ √ N = ⎛⎝ − √ − x x ⎞⎠ − / (4.1)Interestingly, the critical value for x → g c / J = √ N precisely the same as the exact quantumsimulations in . The phase boundary in Eq. (4.1)above is derived as follows: As we see in Fig. 5, forcoupling strengths in the range 1 / < J / g ′ < /√ x m ( J / g ′ ) =[( ( J / g ′ ) − )/ ( J / g ′ ) ] / . This gives rise to a barrierfor the bound states confined within x m ≤ x ≤
1. There-fore, the initial condition determines whether the particleis localized or not by being greater or smaller than x m ; −2 0 2050100 ǫ /g ′ DO S −2 0 2050100 ǫ /g ′ DO S −2 0 2050100 ǫ /g ′ DO S −100 0 100050100 ǫ /g ′ DO S (c)(a) (b)(d) FIG. 7. (Color online) Density of states for different valuesof coupling g . From (a) to (d) J / g ′ = , , . , and 0 . hence, x = x m ( J / g ′ ) gives the expression for the bound-ary between localized and delocalized states in Eq. (4.1).Moreover, when J / g ′ ≤ / x = g ′ √ < (cid:15) < − J + g ′ is classically localized.Let us compare the above result with the quantumspectrum. At the quantum level, we can study the local-ization of eigenstates after adding an infinitesimal sym-metry breaking term H imb = ε ( n L − n R ) where n L ( R ) is polariton number on the left (right) site.In order to compare the quantum results with their clas-sical counterparts, we define an equivalent classical tra-jectory for a given eigenstate of energy (cid:15) n with initialposition (cid:15) n = V l ( x n ) . A comparison with the classicalphase boundary is illustrated in Fig. 6(bottom). For anyvalue of J / g ′ , each eigenstate is represented by a pointwhose horizontal position is x n and whose color reflectsthe quantum expectation value of imbalance defined inEq. (2.1). In this figure, we draw the classical phaseboundary using the criterion whether x n is greater orsmaller than x m . The agreement between the quantumresults and classical picture in polariton bands is quiteremarkable.A signature of the delocalization-localization transi-tion may also be illustrated as changes in the densityof states (DOS), (see Fig. 7). In the weakly interactinglimit g ≪ J , the system is linear and DOS is completelyuniform (Fig. 7(a)). As the coupling g is increased, an-harmonicities emerge in the DOS and eigenvalues startto accumulate near two critical levels (cid:15) c, = g ′ − J and (cid:15) c, = − g ′ + J (Fig. 7(b)). As we continue increasing theratio g / J (Fig. 7(c)), more eigenvalues are depleted attwo regions near touching points of the first and fourthbands with the middle bands, where (cid:15) = ± g ′ √
2, andeventually a gap opens in the DOS (Fig. 7(d)). Usingour WKB quantization conditions near the critical lev-els Eq. (3.4), it is easy to see that the DOS in this region isindeed logarithmically diverging in the system size (totalpolariton number). The level spacing can be genericallycomputed from the quantization condition through δ(cid:15) = (cid:15) n − (cid:15) n − ≃ πN ( d ∆ Sd(cid:15) ) − . Away from the critical levels, the spectrum is harmonicand the density of states is uniform D ( (cid:15) ) ∼ / δ(cid:15) ∼ N .Near the critical points (cid:15) = (cid:15) c,i + λ / N , we have d ∆ S / d(cid:15) ∼ S ( ) ∼ log N (see Appendix E for details) and the DOSscales as D ( (cid:15) ) ∼ N log N that is greater than the regularDOS by a factor of log N . The fact that δ(cid:15) ∼ / N log N also implies that the quantum break times, the time scaleat which the quantum and classical dynamics start todiffer , scale logarithmically with N . This result is alsoconsistent with the N -th order degenerate perturbationtheory for the energy splitting, discussed in , which givesa long-time frequency of ∆ ∼ J ( J / g ) N − for quantum os-cillations, the characteristic time of which scales logarith-mically with N .The anharmoncity observed in the DOS is translatedinto a chaotic behavior at the classical level (see Fig. S3in Appendix B). At the quantum level, the signaturesof quantum chaos can be characterized further in termsof level statistics and the Brody parameter . Thisquantum-to-classical correspondence of the JCDM willbe presented elsewhere. A recent work studies this be-havior for a BEC system where one must consider at leastthree sites (trimer) to make the system non-integrable atthe classical level. Remarkably for JC systems, a two-site (dimer) model would be sufficient for classical non-integrability due to extra degrees of freedom added byqubits. V. DISCUSSION
In summary, we have developed a new technique tostudy the delocalization-localization transition in theJaynes-Cummings dimer model. In this method, weintroduced an effective one-dimensional tight-bindingmodel (polariton basis) for this system and treat eachsite exactly while we treated the intersite polariton hop-ping semiclassically. In order to construct the polari-ton basis, we wrote the Fock space representation forthe JCDM where the Schr¨odinger equation is mappedonto a discrete equation. The continuum approximationcan be made for large system sizes (total polariton num-bers) and the JCDM can be viewed as a moving particleconfined in a potential well details of which are deter-mined by the cavity-qubit interactions. The localizationin this picture is equivalent to the existence of localizedwave-functions near the band edges. We also presented another view for the localization transition in terms ofa gap opening in the density of states. We have foundthat the critical coupling for the transition is the same inboth quantum mechanical and fully classical calculations,which resolves the issue raised by Ref. . Furthermore,we employed a WKB approximation in this picture to de-rive a set of Bohr-Sommerfeld quantization conditions forthe entire energy spectrum in the strong coupling regime.We showed that the quantization rules match with theeigenvalues of the Hamiltonian.The presented method is quite general and can be ap-plied to a variety of models including interactions be-tween spins and bosonic fields. The usual classical limitof these models in the large spin and large photon num-ber limit could involve some difficulties if one wants to in-clude higher order fluctuations due to nonlinear terms fora spin path integral (see Appendix A). Remarkably, ourmethod gets around this difficulty by expanding arounda different classical limit, namely polariton bands. Thepolariton band picture remains valid even at small pho-ton numbers. In particular, this method can be general-ized to other interesting systems where the exact treat-ment of spin-1 / has developed a similar semi-classical approach to study the dynamical transitions dueto quench dynamics in various models including the Dickemodel.Our focus in this article has been the dynamical prop-erties of the JCDM as a closed system. It is worth notingthat the experimental setup is an open system wherephotons can escape from the cavities and the qubits mayrelax and decohere over time. The experimental findingfor the critical coupling is however different from the the-oretical treatments for the closed system. This differencestill remains an open question. In order to accommodatethe dissipation due to environment, it is straightforwardto generalize the classical limit of the JCDM for opensystems using various non-equilibrium path-integral for-malisms originally developed by Kadanoff and Baym or Keldysh . Along these lines, Ref. has recentlystudied the dissipative JCDM by mapping the dynamicsonto a set of Fokker-Planck equations in terms of photonand spin coherent-states in the positive- P representation.They have shown that this method can qualitatively re-produce the experimental measurements. They also con-sidered a driven-dissipative JCDM where it was arguedthat the suppression of the steady-state tunneling currentis a manifestation of the localization transition. It is im-portant to note that the calculations of Ref. have beendone for spin coherent states, which can only be justifiedfor large spins. As we have seen here, treating qubits asspin-1/2 is quite crucial and only this way can one obtainthe full picture. Hence, one interesting future directionis how to generalize the polariton bands to capture thedynamics of an open JCDM.0 ACKNOWLEDGMENTS
This work was initiated by Darius Sadri and has bene-fited immensely from his ideas. Darius tragically passedaway before the work was completed. The authors wouldlike to acknowledge insightful discussions with A. A.Houck, Y. Liu, S. Mandt, J. Raftery, W. E. Shanks, S. J.Srinivasan, H. E. T¨ureci, and D. L. Underwood. In par-ticular, the authors would like to thank D. A. Huse andM. Schir´o for valuable discussions and helpful commentson the manuscript.
Appendix A: Classical dynamics and the localizationtransition
In this appendix, we use the standard coherent-statepath integral formulation of quantum mechanics andwrite down a classical action for the JCDM. We de-rive the classical equations of motion and obtain thelocalization transition after making an analogy with adriven pendulum. In Appendix B, we visualize the time-evolution in the phase space using the Poincar´e sectionsand show that the nonlinear dynamics governed by theclassical equations leads to chaos. The semiclassical pic-ture derived in is based on fully factorizing expectationvalues of products of the photon and the qubit operatorsin Heisenberg equations of motion. Although this typeof factorization is quite customary in the study of JC-based models, there is no physical intuition to what theunderlying assumptions are and it is not clear how onecan systematically improve this process by including newterms.The dynamics of the JCDM with spin- S can be de-scribed in terms of the following action in the real-timecoherent-state path integral formalism S = ∑ s = L,R S sJC [ n s , ¯ ψ s , ψ s ] + S t [ ψ, ¯ ψ ] where the tunneling term becomes S t [ ψ L / R , ¯ ψ L / R ] = J ∫ dt ( ¯ ψ R ψ L + ¯ ψ L ψ R ) (A1)and the JC action is S JC [ n , ¯ ψ, ψ ] =S B [ n , ¯ ψ, ψ ] − ∫ dt ( S n ⋅ B + ν c ¯ ψψ ) , (A2)in which B ( t ) = ν q ˆ z + ig ( ψ − ¯ ψ ) ˆ y + g ( ψ + ¯ ψ ) ˆ x , (A3)here ¯ ψ and ψ are complex fields representing the pho-ton fields and the spin coherent-state is defined in termsof the Bloch state ∣ n ⟩ defined in the coordinate system ( ˆ x , ˆ y , ˆ z ) such that ⟨ n ∣ ˆ S ∣ n ⟩ = S n FIG. S1. The Bloch representation of spin coherent-state. where n is a unit vector (see Fig. S1) spanning the Blochsphere and ˆ S on the LHS is the spin- S operator while S on the RHS denotes the number S . For qubits, we have S = /
2. The first term in the action, Eq. (A2), is theBerry phase contribution S B [ n , ¯ ψ, ψ ] = ∫ dt [ ¯ ψ∂ t ψ + S ⟨ n ∣ ∂ t ∣ n ⟩] . It is important to note that the spin part of the aboveaction is always imaginary even in the imaginary-timeformalism and this signals the fact that this term is topo-logical. Indeed, this term is equal to the area on the Blochsphere enclosed by the path traveled by the Bloch vectorand can be written as ⟨ n ∣ ∂ t ∣ n ⟩ = ∫ dt ∫ dτ n ( t, τ ) ⋅ ( ∂ t n ( t, τ ) × ∂ τ n ( t, τ )) (A4)at which τ is a parameter to define the area on the Blochsphere such that n ( t, ) = n ( t ) n ( t, ) = n here n = ˆ z is the spin quantization direction. In thelimit N → ∞ and S → ∞ , dynamics is described by thestationary solutions of the action, ∂ t n s = B s ( t ) × n s for the qubit and, i∂ t ψ s = ν c ψ s + gS [ n s,x − in s,y ] − Jψ ¯ s for the photon field. Note that the Heisenberg equationsof motion for JCDM Hamiltonian in Eq. (1.2) are ∂ t ˆ S s,k = i [ H, ˆ S s,k ] = ( B s × ˆ S s ) k i∂ t ˆ a s = i [ H, ˆ a s ] = ν c ˆ a s + g [ ˆ S x − i ˆ S y ] where B is defined in Eq. (A3) and the additional sub-script k denotes the k -th component of a vector. Theclassical equations of motion are equivalent to factoriz-ing approximation of the expectation values of quantumoperators; i.e. ∂ t ⟨ ˆ S ⟩ ≈ ⟨ B ⟩ × ⟨ ˆ S ⟩ ⟨ ... ⟩ denotes the expectation value over the initialstate and we make the following identifications: ⟨ a ⟩ = ψ and ⟨ ˆ S ⟩ = S n .It is more convenient to represent spin in the sphericalcoordinate system n = ( sin θ cos φ, sin θ sin φ, − cos θ ) , asin Fig. S1, and to split the cavity fields into their real andimaginary parts ψ s = R s + iI s . Using this parametriza-tion, the equations of motion become˙ φ s = − ν q − g ( R s cos φ s − I s sin φ s ) cot θ s (A5a)˙ θ s = g ( R s sin φ s + I s cos φ s ) (A5b)for the spin and˙ R s = ν c I s − gS sin θ s sin φ s − JI ¯ s (A5c)˙ I s = − ν c R s − gS sin θ s cos φ s + JR ¯ s (A5d)for the cavity fields, where s = L, R and its opposite ¯ s = R, L . Here, our focus is the resonant case where ν q = ν c and the free dynamics of the spin and photon field canbe removed in the rotating frame.The classical dynamics of the JCDM takes place inan eight dimensional phase space. In order to finda lower bound on the critical value of g / J where thedelocalization-localization transition occurs, we approachthe critical point from the localized side. In the localizedphase, we start with all polaritons on the left cavity andwe can assume ψ L ( t ) ≈ √ N and ψ R ( t ) ≈ N ≫ S . So, thedynamics in the rotating frame simplifies into˙ R L = − gS sin θ L sin φ L , (A6a)˙ R R = − gS sin θ R sin φ R , (A6b)˙ I L = − gS sin θ L cos φ L , (A6c)˙ I R = − gS sin θ R cos φ R + J √ N , (A6d)where only the fourth equation couples the two JC is-lands. This equation describes a motion subject toan external drive F = J √ N . Depending on the ra-tio η = gS / J √ N , the sign of the RHS may be allowedto change over time or not. These two modes of mo-tion correspond to delocalized and localized phases. If η <
1, the RHS of Eq. (A6d) always remains positiveand the right (empty) cavity absorbs more and more en-ergy. Therefore, excitations can be fully transferred tothe empty cavity and this process continues repeatedly.However, if there exists a turning point at which the RHSof Eq. (A6d) vanishes (i.e. η sin θ cos φ = η c =
1. So, the lowest value of the θ R ( t = 0) / π g ( c l ) c / J √ N FIG. S2. (Color online) The classical critical ratio g / J asa function of initial configuration for the spin. Here, the dy-namics is given by Eqs. (A5a-d) and the other spin anglesare initially fixed at θ L ( ) = φ L ( ) =
0. The green curve isgiven by Eq. (A9), based on analogy to the pendulum, andthe blue circles is the critical boundary obtained numericallyafter calculating the long-time averaged imbalance. We put S = / coupling ratio g / J to get a localized phase is given by g ( cl ) c J = S √ N . (A7)After plugging S = /
2, we see that this result is the sameas the quantum simulations as well as the WKB pic-ture studied in this article. However, this is differentfrom the previous numerical analysis on the semiclassicalfactorized Heisenberg equations of motion where thecritical ratio is found to be g / J ≈ √ N that is off bya factor of √ only considered the initial con-dition with both qubits down and ignored the fact thatthe critical ratio g / J actually depends on the initial statesof qubits. Figure S2 illustrates this dependence explicitly.As we see in this plot, the phase boundary in terms of g / J √ N starts around 2 √ ≃ .
8, the same as Ref. , butreaches a minimum at 2 for θ R = π /
2. Hence, the classi-cal dynamics (i.e. the dynamics in the limit N ≫ S ≫ g required to obtainself-trapping is given by Eq. (A7).Few remarks regarding Fig. S2 are in order. First,note that our analysis of the critical coupling for arbi-trary initial spin configurations is based on the dynamicsdescribed by Eqs. (A6a-d) where the initial state of theleft qubit θ L and φ L corresponding to the fully populatedJC site does not make any significant change and henceis fixed at θ L ( ) = φ L ( ) =
0. Second, since we are look-ing for a lower bound on the critical ratio g / J , we choose φ R ( ) = g in thefirst term at the RHS of Eq. (A6d). Given these initialconditions, the dynamics of I R and θ R can be combined2into a single equation˙ θ R = g S ( cos θ R − cos θ R ( )) + Jg √ N ( θ R − θ R ( )) (A8)where the initial conditions are θ R ( ) ≠ θ R ( ) = θ R ) and subjectto a constant torque τ = Jg √ N . It is easier to un-derstand the localization transition using the pendulumanalogy. The pendulum may do a full rotation (delocal-ized) or only oscillate with a small amplitude (localized).In the oscillatory mode the motion is bound by the turn-ing points at which ˙ θ R =
0. So, the existence of turningpoints can be used to determine the boundary betweenlocalized and delocalized dynamics. The critical valuecan then be written as g ( cl ) c / J √ N = / sin θ where θ isthe non-zero solution ofsin θ + cos θ − cos θ R ( ) θ − θ R ( ) = . (A9)Notice that the critical ratio for θ R ( ) ≥ π / g ( cl ) c / J √ N = / sin θ R ( ) as expected from settingthe RHS of Eq. (A6d) zero. However, for θ R ( ) < π / g ( cl ) c / J √ N = / sin θ R ( ) . This means thatthe spin dynamics must be considered fully along withthe cavity field dynamics to yield the correct criticalvalue. Finally, we confirm our analysis by numericallysimulating the full nonlinear classical dynamics (shownas blue circles in Fig. S2) and taking the averaged im-balance as an order parameter to determine the phaseboundary between localized and delocalized phases. Notethat in numerics we only set the initial condition for theleft spin θ L ( ) = φ L ( ) = q i ( i = q i ( t ) = q ( cl ) i ( t ) + r i ( t ) can be approximated by S[ q ] ≈ S[ q ( cl ) ] + ∫ dt ′ ∫ dt ′′ r i ( t ′ ) G ij ( t ′ , t ′′ ) r j ( t ′′ ) + ⋯ where the classical solutions satisfy the stationary con-dition δ S/ δq ( cl ) i = G ij ( t ′ , t ′′ ) = δ S[ q ] δq i ( t ′ ) δq j ( t ′′ ) ∣ q i = q ( cl ) i . The difficulty in calculating the Green’s function is dueto the Berry phase term for spins in Eq. (A4) which is notquadratic and causes G ij ( t ′ , t ′′ ) to explicitly depend on the time variables t ′ and t ′′ during a trajectory. There-fore, finding the Green’s function for an arbitrary trajec-tory is generically a difficult task especially since there isno close solution to Eqs. (A5(a)-(d)). This process canhowever be done numerically.The coupling to the qubit not only induces a non-linearity that leads to the localization, but also enlargesthe phase space available to the system by adding moredegrees of freedom. Unlike the BECDW , the JCDM isnot integrable for non-zero coupling. The study of chaosin the JCDM is beyond the scope of the present paper. Appendix B: Classical chaos and the KAM theorem
As mentioned in Appendix A, the non-linearity due tocoupling to the qubit may lead to chaos at the classicallevel. When g =
0, our system is (Liouville) integrable ,periodic or quasi-periodic motion then occurs on invari-ant tori with angles as variables and different tori labeledby conserved action I i = π ∮ γ i p ⋅ dq . (B1)The Kolmogorov-Arnold-Moser (KAM) theorem statesthat as a non-linear perturbation is introduced, the ratio-nal (resonant/periodic) tori are destroyed and the quasi-periodic ones are deformed. As the interaction is in-creased, only sufficiently irrational tori continue to sur-vive, those that do form a Cantor set. Let us now visual-ize this behavior in our system. We use the restricted dy-namics to construct Poincar´e sections. For certain initialconditions, the dynamics is restricted to a four dimen-sional subspace. In this construction a 4d phase spacegives rise to a 3d fixed energy surface. To construct thePoincar´e section, one follows the dynamics and recordsthe state of the system when one of the degrees of free-dom, called periodic variable, reaches a certain value.This specifies a point in the 2d plane corresponding tothe other degrees of freedom. The section is then filledin by sampling initial conditions on the fixed energy sur-face and following the dynamics for each initial state fora long time.An example of such sections is given in Fig. S3 for var-ious values of g / J . Here, the initial conditions are takenas I L = R R = φ L = π /
2, and φ R = Z ≡ n L − n R = N for R L = √ N and I R = t = θ L = gR L (B2a)˙ θ R = gI R (B2b)˙ R L = − gS sin θ L − JI R (B2c)˙ I R = − gS sin θ R + JR L (B2d)3 FIG. S3. (Color online) Poincar´e sections for various val-ues of g / J √ N . (a)-(d) 3D diagrams showing explicitlyone of spin angles, (e)-(h) 2D diagrams (equivalent to thered sections of its left figure) using the polar coordinates ( r, α ) = (√( R L + I R )/ N, tan − ( I R / R L )) . From top to bot-tom g / J √ N = . , . , . , and 1 . together with ˙ R R = ˙ I L = ˙ φ L = ˙ φ R = . (B3)In this four dimensional subspace the tunneling potentialin Eq. (A1) is identically zero, which means the energyassociated with this subspace is independent of the hop-ping amplitude J . In Fig. S3, the periodic variable ischosen to be θ L ( t ) (with 2 π period) and the other twovariables are the cavity fields R L ( t ) and I R ( t ) . As wesee, for sufficiently small values of g / J as in Fig. S3 (a)and (e) the motion is periodic while it becomes more andmore chaotic as the ratio g / J is increased. Nevertheless,there are few surviving quasi-periodic orbits for interme-diate interactions.Interestingly, the restricted equations of motion westudy for classical chaos are identical to the classicalequations of motion for electrons moving in a 2d peri-odic potential subject to an external magnetic field. Thissystem has been previously studied and is known todisplay a variety of ballistic, 1d and 2d normal diffusiveand anomalous diffusive transport, related to the chaoticbehavior. Appendix C: Review of the Bose-Einsteincondensate double-well problem
The BECDW system is considerably simpler than theJCDM as there is no spin degrees of freedom. TheSchr¨odinger equation has the same form as in Eq. (2.6)with two differences that the hopping term is B ( x ) = J / √( + x + h )( − x + h ) and the nonlinearity is ex-plicit in the potential energy term W ( x ) = γx + εx where γ is the interaction strength and ε is an ad hoc term to break the left-right parity symmetry. We choosethe ansatz Ψ ( x, t ) = α ( x, t ) e ih S ( x,t ) with a real-valuedphase and amplitude. The real and imaginary parts ofthe Schr¨odinger equation up to first order in h are foundto be − ∂ t S = − B cos ( ∂ x S ) + W ( x ) − hB cos ( ∂ x S )− ∂ t α = αB cos ( ∂ x S ) ∂ x S + ( αB ′ + α ′ B ) sin ( ∂ x S ) (C1)where B ( x ) = J √ − x B ( x ) = J √ − x . In order to identify the Hamiltonian and the canonicalmomentum, we use the Hamilton-Jacobi equations H = − ∂ t S ϕ = ∂ x S So, the Hamiltonian can be written as H = − B ( x ) cos ( ϕ ) + W ( x ) + W Q ( x ) quantum potential , W Q ( x, ϕ ) = − hB ( x ) cos ( ϕ ) + O ( h ) . We introduce the classical velocity v = ∂H∂ϕ = B ( x ) sin 2 ϕ, hence, the second identity in Eq. (C1) can be rewrittenin terms of the probability density P = α , as ∂ t P + ∂ x ( P v ) = . which is the continuity equation. Now that the Hamilto-nian is determined, one can solve for the stationary solu-tions Ψ ( x, t ) = α ( x ) e − ih (cid:15)t e ih S ( x ) where S ( x ) is expandedin powers of h , S ( x ) = S ( ) ( x ) + hS ( ) ( x ) + O ( h ) , thus, the canonical momentum can be also expanded as ϕ = ∂ x S = ϕ ( ) + hϕ ( ) + O ( h ) . Next, we can solve for phase space trajectories with afixed energy (cid:15) by plugging this expansion into the Hamil-tonian (cid:15) = − B cos ( ϕ ( ) + hϕ ( ) ) + W ( x ) + W Q ( x, ϕ ( ) ) At zeroth order, we get ϕ ( ) ( x ) =
12 Arccos ( (cid:15) − W ( x )− B ( x ) ) and the first order correction would be ϕ ( ) = B ( x ) B ( x ) cos ( ϕ ( ) ( x )) sin ( ϕ ( ) ( x ))= − ( − x ) (cid:15) − W ( x )( J ( − x ) − ( (cid:15) − W ( x )) ) / . Using these results, quantization rules can be found ac-cordingly . Appendix D: Exact solutions of the Schr¨odingerequation near boundaries
Here, we briefly discuss the exact solutions at theboundaries we encounter in JCDM polariton bands. Ex-cept for the critical region which will be explained at theend, the linearized Schr¨odinger equation around the clas-sical turning point x = z α + ξ up to O ( h ) corrections canbe easily shown to be ∂ ξ ψ i − β α ∂ ξ ψ i − c i,α γ α ξ ψ i ≈ , (D1) where z α is the turning point that is the solution to (cid:15) = V αi ( z α ) , in which α = l, h refers to the stationary pointswith ϕ i = , π / β α = z α − z α ,γ α = √ − z α h J ∣ ∂ x V αi ( z α )∣ , here, a sign coefficient c i,α is introduced to be ± ξ > ∣ z α ∣ = J (( (cid:15) + J − g ′ ) J − ( g ′ ∓ (cid:15)J ) + ∣ g ′ ∓ (cid:15)J ∣√ g ′ + J ∓ (cid:15)J ) / . (D2)where the ± signs are both accepted only if there arefour turning points on both inner and outer boundaries.The solution of the linearized differential equation can bewritten analytically in terms of the Airy functions ψ i = e β α ξ / [ C a Ai ( c i,α γ / α ˜ ξ ) + C b Bi ( c i,α γ / α ˜ ξ )] where ˜ ξ = ξ + c i,α β α / γ α . It is important to note thatin deriving the above equation near the turning pointat which ϕ i → π /
2, the wave-function is oscillating sofast ψ i ∝ e iπx / h ; hence, we write down the linearizedSchr¨odinger equation for the slowly varying part of theactual wave-function by mutliplying the wave-functionby e − iπx / h to cancel the oscillating factors.In the case of upper-upper or lower-lower bands, thereis a critical region, near energy (cid:15) c, = g ′ − J or (cid:15) c, =− g ′ + J respectively, between the localized and delocal-ized states where the potential term V l ( x ) or V h ( x ) isextremum at the turning point x =
0. This situationrequires a separate analysis as the exact solution is nolonger the Airy functions and the quantization conditionwill be different from the usual Bohr-Sommerfeld quan-tization. Near the critical level, energy can be written as (cid:15) = (cid:15) c,k + λ / N where the expanded Schr¨odinger equationnear x = + ξ is simplified into ∂ ξ ψ k − ξ ∂ ξ ψ k + c k Jh [− λh + ξ ( g ′ − J )] ψ k ≈ c k is c = − c = + , ψ k = e ξ / [ C a D iχ − / ( c k e iπ / ξ √ µ / h )+ C b D iχ − / ( c k e − i π / ξ √ µ / h )] , where we neglect the terms of order unity compared to1 / h and µ = √ g ′ J − , χ = λ µJ . (D4)5 Appendix E: Derivation of quantization conditionsfor the lower-lower band
Here we explicitly show how to derive the quantizationrules for lower-lower polaritons. (i) Delocalized states, (cid:15) < J − g ′ :Here, the particle is confined in a two-sided poten-tial. Since the classically allowed region is defined be-tween two turning points ∣ x ∣ < z l , the WKB wave-function(Fig. 3(b)) would simply be ψ ( x ) = ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ α ( x ) ( C + e ih S ( x ) + C − e − ih S ( x ) ) , ∣ x ∣ < z l α ( x ) C ′− e − h Q ( x ) , x > z l α ( x ) C ′+ e h Q ( x ) . x < − z l (E1)and the other modes are described by ψ ( x ) = α ( x ) e − h Q ( x ) ,ψ ( x ) = α ( x ) e h Q ( x ) . The above modes will add a correction of order h to thewave-function exponents in Eq. (E1) through Eq. (2.14);however, these terms are exponentially small and will beneglected. Matching the WKB wave-function with theAiry functions at x = ± z l leads to the following connectionformulas i C − C + = e ih ∆ S , C b = S = ∫ z l ϕ dx . Similarly, one can derive theconnection formula for x = − z l + ξi C + C − = e ih ∆ S . Thus, we get h ∆ S = ( n + ) π which means1 h ∫ z l − z l ϕ ( x ) dx = ( n + ) π (E2)This is similar to the familiar Bohr-Sommerfeld quanti-zation for a particle confined in a potential well . (ii) Localized states, J − g ′ < (cid:15) < − g ′ √ x = ±
1. Hence, there are two setsof turning points z l and z h : One is bouncing off the ϕ = (cid:15) = V l ( x = z l ) , and the other is bouncing off the 2 ϕ = π potential,which is the solution of (cid:15) = V h ( x = z h ) . Provided thatthe classically allowed region is z h < ∣ x ∣ < z l , the WKBwave-function (Fig. 3(c)) would be ψ ( x ) = ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ α ( x ) C ′− e − h Q ( x ) , x > z l α ( x ) ( C R, + e ih S ( x ) + C R, − e − ih S ( x ) ) , z h < x < z l α ( x ) ( C ′′+ e h Q ( x ) + C ′′− e − h Q ( x ) ) , ∣ x ∣ < z h α ( x ) ( C L, + e ih S ( x ) + C L, − e − ih S ( x ) ) , − z l < x < − z h α ( x ) C ′+ e h Q ( x ) . x < − z l Close to ∣ x ∣ = z l the connection formula is the same asthe previous case i C R, − C R, + = e ih ∆ S ,i C L, + C L, − = e ih ∆ S (E3)where the boundary of the integral is only slightly differ-ent as in ∆ S = ∫ z l y r ϕ dx and y r is the solution to theequation (cid:15) = W ( x, ϕ ) corresponding to ϕ = π / y r = ∣ (cid:15) ∣ g ′ √ − ( (cid:15) g ′ ) . (E4)Near x = z h , ϕ → π and the wave-function is oscillat-ing so fast and ψ ∝ e iπx / h ; therefore, as mentioned inAppendix D, only the slowly varying part of the wave-function will be matched to the Airy functions and therewill be an additional term πy r / h in the quantization ruleas a result of this manipulation.Close to x = z h and on the forbidden side, the connec-tion formula yields C a C b = e h ∆ Q (E5)where ∆ Q = ∫ z h − z h φ dx . The matching condition for theallowed side is simplified into RRRRRRRRRRR C R, + e ih ∆ S ′ − iπy r / h + iπ / + C R, − e − ih ∆ S ′ + iπy r / h − iπ / − C R, + e ih ∆ S ′ − iπy r / h + iπ / + C R, − e − ih ∆ S ′ + iπy r / h − iπ / RRRRRRRRRRR = ∣ C b C a ∣ , where ∆ S ′ = ∫ y r z h ˜ ϕ dx , and˜ ϕ ( x ) =
12 Arccos ⎛⎝ (cid:15) − W ( x ) B ( ) ( x ′ ) ⎞⎠ . (E6)It is worth noting that there is no minus sign inside theargument as opposed to Eq. (2.13). Combining this withEq. (E3), we get1 − cos 2∆ S (eff)4 + cos 2∆ S (eff)4 = e − h ∆ Q where ∆ S (eff)4 = ∆ S − ∆ S ′ + πy r / h = ∫ z l y r ϕ dx −∫ y r z h ˜ ϕ dx + πy r / h . This can be recast in the followingquantization condition1 h ∆ S (eff)4 ( (cid:15) ) = nπ ± e − h ∆ Q ( (cid:15) ) . (E7)As a result, the tunneling rate is related to the energysplitting in the energy spectrum through δ(cid:15) (split) = e − h ∆ Q ( (cid:15) ) RRRRRRRRRRR d ∆ S (eff)4 d(cid:15) RRRRRRRRRRR − . (iii) Critical region: (cid:15) ∼ J − g ′ :Consider a small deviation from the critical energy (cid:15) = J − g ′ + λ / N . The WKB solutions close to the criticalpoint x = ∂ x S ( ) = π − µx ,∂ x S ( ) = λ Jµ x ,α ( x ) = α , √ µJx e x / . (E8)So, the wave-function becomes ψ ( x ) = C ± exp (± iπ h x ∓ i µx h + (± iχ − / ) log ( x √ µ / h )) (E9)and the quantization formula is derived to bearg [ e − ih ∆ S ( eff ) √ π Γ ( / − iχ ) − e πχ / ] = nπ + π / S ( eff ) = ∫ z l y r ϕ dx − ∫ y r ˜ ϕ dx + πy r , (E11) χ and y r are defined in Eqs. (D4) and (E4), respectively. Appendix F: Quantization conditions for otherpolariton bands
In this appendix, we show the quantization conditionsfor upper-upper and middle polariton bands.
Upper-upper band:
The upper-upper mode is pop-ulated for energies g ′ √ ≤ (cid:15) ≤ J + g ′ , as illustrated inFig. 2. The situation is very similar to the lower-lowerpolariton band with the only difference that the localizedstates appear in the bottom of the band and the delocal-ized states appear at higher energies. Likewise, there arethree regions: one with delocalized wave-functions cen-tering around x =
0, one with localized wave-functionsin which the probability density is maximized close tothe boundary points x = ±
1, and one at the critical re-gion separating the localized and delocalized states. Thisband is similar to the BECDW with attractive interac-tion. The derivation of the quantization rules is verysimilar to Appendix E and we shall only quote the finalresults. (i) Delocalized states, (cid:15) > − J + g ′ : The particle is con-fined in a two-sided potential well and the wave-functionis large close to the center ( x = ϕ = π and the boundary match-ing must be done for a slowly varying component of thewave-function. The quantization condition is found by1 h ∫ z l − z l ˜ ϕ ( x ) dx = ( n + ) π (F1) where ˜ ϕ ( x ) =
12 Arccos ⎛⎝ (cid:15) − W ( x ) B ( ) ( x ′ ) ⎞⎠ . (F2)Notice that there is no minus sign in the denominator asopposed to Eq. (2.13) which has to do with only takingthe slowly varying multiple of the wave-function. (ii) Localized states, g ′ √ < (cid:15) < − J + g ′ : The potentiallooks like a double-well with two minima at x = ±
1. Thequantization rule would be1 h ∆ S (eff)1 ( (cid:15) ) = nπ ± e − h ∆ Q ( (cid:15) ) (F3)where ∆ S (eff)1 = ∫ z l y r ˜ ϕ dx − ∫ y r z h ϕ dx + π y r , ∆ Q = ∫ z h − z h φ dx, in which ˜ ϕ is defined in Eq. (F2) and other parametersare introduced as in Appendix E. (iii) Critical region, (cid:15) ∼ − J + g ′ : Without much dif-ference from the lower-lower band, here the quantizationformula is found to bearg [ e − ih ∆ S ( eff ) √ π Γ ( / + iχ ) − e − πχ / ] = nπ + π / S ( eff ) = ∫ z l y r ˜ ϕ dx − ∫ y r ϕ dx + πy r ,χ and y r are defined in Eqs. (D4) and (E4), respectively. Middle bands:
The middle polariton bands cor-respond to the energy range ∣ (cid:15) ∣ ≤ g ′ √
2. The wave-function has two non-zero components: the lower-upperand upper-lower polariton modes. This region in turn canbe divided to two subregions: In the first region, when J ≤ (cid:15) ≤ g ′ √ − g ′ √ ≤ (cid:15) ≤ − J , one component is domi-nant in each half ( x > x <
0) since only one compo-nent is classically allowed and the other one is forbidden.In the second region, when ∣ (cid:15) ∣ < J both components arenon-zero at the same region of x . However, in terms ofquantization relations the treatments for both cases areidentical as the first order correction, in Eq. (2.10), doesnot couple the second and third components directly andall corrections come in as higher order contributions (atleast of order h ). Therefore, it does not really matterwhether these two components overlap over a range of x or not.The upper-lower (second) polariton component is con-fined between two turning points z l ≤ x ≤ z h where z l and z h are solutions to (cid:15) = V l ( x ) and (cid:15) = V h ( x ) givenin Eq. (D2). The lower-upper (third) polariton compo-nent is then confined in − z h ≤ x ≤ − z l due to symmetry7(see Fig. 2). So, the WKB wave-function takes the formΨ = [ , ψ , ψ , ] T where ψ ( x ) = ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ α ( x ) ( C , + e ih S ( x ) + C , − e − ih S ( x ) ) , z l < x < z h α ( x ) C ′ , − e − h Q ( x ) , x > z h α ( x ) C ′ , + e h Q ( x ) , x < z l and ψ ( x ) = ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ α ( x ) ( C , + e ih S ( x ) + C , − e − ih S ( x ) ) , − z h < x < − z l α ( x ) C ′ , − e − h Q ( x ) , x > − z l α ( x ) C ′ , + e h Q ( x ) . x < − z h Starting from the connection formulas at the two bound-aries for each component, it is straightforward to derivethe quantization conditions1 h ∆ S (eff)2 ( (cid:15) ) = nπ − π y r (F5a)1 h ∆ S (eff)3 ( (cid:15) ) = nπ + π y r (F5b)where ∆ S (eff)2 = ∫ z l y r ϕ dx − ∫ y r z h ˜ ϕ dx ∆ S (eff)3 = ∫ − z l − y r ϕ dx − ∫ − z h − y r ˜ ϕ dx and y r is defined in Eq. (E4). It is easy to show thatthese two relations are indeed identical as the integrandsare even functions of x . Appendix G: Validity of WKB analysis
The derivation of WKB relies on our original assump-tion that the total polariton number N is large and a1 / N expansion is applicable. A natural question is towhat extent this assumption can be justified. In order tocheck the validity of the WKB approximation, we use theHusimi-Kano Q representation to compare our semi-classical picture and the quantum eigenstates (see a com-prehensive discussion of the Husimi function in Ref. ).Similar analyses have been done for the BEC double-wellsystem .As we have seen in Sec. II, dynamics within i -th polari-ton band can be described by two conjugate variables the“position” x = Z / N and the “momentum” ϕ i . Let us fo-cus on the lower-lower polariton band and define θ = ϕ .The Poisson bracket is promoted to the canonical com-mutation relation [ x, θ ] = ih at the quantum level wherethe Planck’s constant is h = / N . A squeezed coherent-state can be represented in position and momentum, ⟨ θ ′ ∣ θ + ix ⟩ = ( πκ ) / exp (− i xθ ′ h − ( θ ′ − θ ) κ )⟨ x ′ ∣ θ + ix ⟩ = ( π / κ ) / exp ( i x ′ θh − κ ( x ′ − x ) ) . (G1) x θ / π −1 −0.5 0 0.5 1−1−0.500.51(a) FIG. S4. (Color online) Comparison of the classical phasespace and the quantum Husimi distribution in the lower-lowerpolariton band. (a) classical energy contours. The Q distri-bution is shown (b) for the ground state, (c) for the mid band(oscillatory) state 52-th, (d) close to the critical level (separa-trix in (a)) 73-rd, (e) for the localized state 97-th state. Here N =
100 and J / g ′ = / The squeezing parameter κ is the key parameter inconnecting the classical picture and quantum eigen-states . It determines the relative resolution in thephase space ( x, θ ) . The optimum κ generally dependson the explicit form of the Hamiltonian in terms of theconjugate variables of interest. For a simple harmonicoscillator, this parameter must be chosen equal to thezero-point fluctuations κ = ( mω /̵ h ) / in ( x, p ) -space andequal to one κ = ( a, a † ) -space (theusual Q function in quantum optics ). For our case, κ must be tuned to an optimum value to obtain the bestresolution. In general, we write κ = sκ where s > s = κ = ( N g / J ) / comes from the harmonicapproximation of the Hamiltonian close to the ground8going up in spectrum d ec r e a s i n g N FIG. S5. (Color online) Husimi Q representation for various values of system size N in the lower-lower polariton band: (a) N =
20 (b) N =
10, and (c) N =
6. From left to right, each row shows the ground state, a mid-band (oscillatory) state, a stateclose to the critical level, and a localized state. Here J / g ′ = / state H ≈ − g ′ − J + J θ + g ′ x . (G2)For a pure state ∣ ψ ⟩ , the Husimi function is defined by Q ( x, θ ) = ∣⟨ θ + ix ∣ ψ ⟩∣ (G3)In order to compute the inner product numerically, weuse the following identity ⟨ θ + ix ∣ ψ ⟩ = ( π / κ ) / ∑ x ′ =− C ( x ′ ) exp ( i x ′ θh − κ ( x ′ − x ) ) (G4)where C ( Z ) is the fourth-component of the wave-function in the polariton basis. Note that in these cal-culations only C ( Z ) is taken into account as the rest ofthe components are negligible in the lower-lower band.Figure S4 shows that the Husimi functions of the Hamil-tonian eigenstates match quite well with the classical phase space energy contours for large values of N . Aswe see in Fig. S4(b), the ground state is represented bya point-like distribution meaning that it is a minimumuncertainty (Gaussian) wave-packet, consistent with ourharmonic approximation above. The next excited states(Fig. S4(c)) are oscillatory states. At the critical level (cid:15) c, = − g ′ + J , we arrive at the separatrix (bifurcationpoint) of the classical phase space, which is also manifestin the Husimi function as in Fig. S4(d). Ultimately, abovethe critical level we obtain the localized states which arerepresented by two seperate branches close to x = − N =
20, 10, and 6 inFigs. S5(a)-(c). The classical phase space is independentof N and is shown in Fig. S4(a). It is clear that the9 ǫ /g ′ -3 -2 -1 0 1 2 3 E rr o r i n W K B FIG. S6. (Color online) Error in the Bohr-Sommerfeld quan-tization rules in the entire spectrum. Different symbols referto different quantization conditions and the vertical lines showvarious energy scales; both are explained in the caption ofFig. 4. Here J / g ′ = N = fluctuations are more pronounced as N becomes smaller;however, even up to such small values as N = J / g ′ is illustrated in Fig. S6. † Deceased. A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. Hol-lenberg, Nat Phys , 856 (2006). M. J. Hartmann, F. G. S. L. Brandao, and M. B. Plenio,Nat Phys , 849 (2006). M. Hartmann, F. Brando, and M. Plenio, Laser and Pho-tonics Reviews , 527 (2008). A. Tomadin and R. Fazio, J. Opt. Soc. Am. B , A130(2010). A. A. Houck, H. E. T¨ureci, and J. Koch, Nat Phys , 292(2012). S. Schmidt and J. Koch, Annalen der Physik , 395(2013). D. G. Angelakis, M. Huo, E. Kyoseva, and L. C. Kwek,Phys. Rev. Lett. , 153601 (2011). M. Schir´o, M. Bordyuh, B. ¨Oztop, and H. E. T¨ureci, Phys.Rev. Lett. , 053601 (2012). J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hartmann,Phys. Rev. Lett. , 163605 (2013). J. Otterbach, M. Moos, D. Muth, and M. Fleischhauer,Phys. Rev. Lett. , 113001 (2013). A. Le Boit´e, G. Orso, and C. Ciuti, Phys. Rev. Lett. ,233601 (2013). F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E. T¨ureci,and J. Keeling, Phys. Rev. Lett. , 233603 (2012). I. Carusotto and C. Ciuti, Rev. Mod. Phys. , 299 (2013). K. Le Hur, L. Henriet, A. Petrescu, K. Plekhanov,G. Roux, and M. Schir´o, arXiv:1505.00167 (2015). J. Koch, A. A. Houck, K. L. Hur, and S. M. Girvin, Phys.Rev. A , 043811 (2010). A. Petrescu, A. A. Houck, and K. Le Hur, Phys. Rev. A , 053804 (2012). J. Cho, D. G. Angelakis, and S. Bose, Phys. Rev. Lett. , 246809 (2008). M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor,Nat Photon , 1001 (2013). R. O. Umucalilar and I. Carusotto, Phys. Rev. A ,043804 (2011). M. Hafezi, M. D. Lukin, and J. M. Taylor, New J. Phys. , 063001 (2013). M. Brune, F. Schmidt-Kaler, A. Maali, J. Dreyer, E. Ha-gley, J. M. Raimond, and S. Haroche, Phys. Rev. Lett. , 1800 (1996). K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E.Northup, and H. J. Kimble, Nature , 87 (2005). J. P. Reithmaier, G. Sek, A. Loffler, C. Hofmann, S. Kuhn,S. Reitzenstein, L. V. Keldysh, V. D. Kulakovskii, T. L.Reinecke, and A. Forchel, Nature , 197 (2004). T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova, H. M.Gibbs, G. Rupper, C. Ell, O. B. Shchekin, and D. G.Deppe, Nature , 200 (2004). E. Peter, P. Senellart, D. Martrou, A. Lemaˆıtre, J. Hours,J. M. G´erard, and J. Bloch, Phys. Rev. Lett. , 067401(2005). A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S.Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J.Schoelkopf, Nature , 162 (2004). A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J.Schoelkopf, Phys. Rev. A , 062320 (2004). S. Schmidt, D. Gerace, A. A. Houck, G. Blatter, and H. E.T¨ureci, Phys. Rev. B , 100507 (2010). G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls,Phys. Rev. A , 4318 (1997). A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy,Phys. Rev. Lett. , 4950 (1997). M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani,and M. K. Oberthaler, Phys. Rev. Lett. , 010402 (2005). S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Na-ture , 579 (2007). J. Esteve, C. Gross, A. Weller, S. Giovanazzi, and M. K.Oberthaler, Nature , 1216 (2008). M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov,H. Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaitre,G. Malpuech, and J. Bloch, Nat Phys , 275 (2013). J. Raftery, D. Sadri, S. Schmidt, H. E. T¨ureci, and A. A.Houck, Phys. Rev. X , 031043 (2014). P. A. Braun, Rev. Mod. Phys. , 115 (1993). B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov,Phys. Rev. Lett. , 2803 (1997). M. J¨a¨askel¨ainen and P. Meystre, Phys. Rev. A , 043603(2005). B. Juli´a-D´ıaz, D. Dagnino, M. Lewenstein, J. Martorell,and A. Polls, Phys. Rev. A , 023615 (2010). L. Fu and J. Liu, Phys. Rev. A , 063614 (2006). E. M. Graefe and H. J. Korsch, Phys. Rev. A , 032116(2007). V. S. Shchesnovich and M. Trippenbach, Phys. Rev. A ,023611 (2008). F. Nissen and J. Keeling, Phys. Rev. A , 063628 (2010). L. E. Ballentine,
Qunatum Mechanics A Modern Develop-ment , 1st ed. (World Scientific Publishing Co., Singapore,1998). T. Brody, Lettere Al Nuovo Cimento Series 2 , 482 (1973). V. Manfredi, Lettere al Nuovo Cimento , 135 (1984). I. Tikhonenkov, A. Vardi, J. R. Anglin, and D. Cohen,Phys. Rev. Lett. , 050401 (2013). B. Sciolla and G. Biroli, Phys. Rev. B , 201110 (2013). L. P. Kadanoff and G. Baym,
Quantum Statistical Mechan-ics (W. A. Benjamin, Inc., New York, 1962). A. Kamenev,
Field Theory of Non-equilibrium Systems (Cambridge University Press, Cambridge, 2011). S. Mandt, D. Sadri, A. A. Houck, and H. E. T¨ureci, NewJ. Phys. , 053018 (2015). E. Fradkin,
Field Theories of Condensed Matter Physics ,2nd ed. (Cambridge University Press, Cambridge, 2013). A. Altland and B. Simons,
Condensed Matter Field Theory ,2nd ed. (Cambridge University Press, Cambridge, 2010). S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy,Phys. Rev. A , 620 (1999). J. V. Jos´e and E. J. Saletan,
Classical dynamics a contem-porary approach (Cambridge University Press, Cambridge,1998). T. Geisel, J. Wagenhuber, P. Niebauer, and G. Obermair,Phys. Rev. Lett. , 1581 (1990). E. Whittaker and G. Watson,
A Course of Modern Anal-ysis (Cambridge University Press, Cambridge, 1927). K. Husimi, Proceedings of the Physico-Mathematical So-ciety of Japan. 3rd Series , 264 (1940). H.-W. Lee, Physics Reports , 147 (1995). K. W. Mahmud, H. Perry, and W. P. Reinhardt, Phys.Rev. A , 023615 (2005). M. O. Scully and M. S. Zubairy,