Synthetic gauge fields and homodyne transmission in Jaynes-Cummings lattices

Many-body physics is traditionally concerned with systems of interacting massive particles. Recent studies of effective interactions between photons, induced in the circuit QED architecture by coupling the microwave field to superconducting qubits, have paved the way for photon-based many-body physics. We derive the magnitude and intrinsic signs of photon hopping amplitudes in such circuit QED arrays. For a finite, ring-shaped Jaynes-Cummings lattice exposed to a synthetic gauge field we show that degeneracies in the single-excitation spectrum emerge, which can give rise to strong correlations for the interacting system with multiple excitations. We calculate the homodyne transmission for such a device, explain the generalization of vacuum Rabi splittings known for the single-site Jaynes-Cummings model, and identify fingerprints of interactions beyond the linear response regime.


Introduction
The notion of a quantum simulator can be traced back to a remark by Richard Feynman [1] which was later made more precise by David Deutsch [2] and Seth Lloyd [3]. Since simulating quantum many-body systems on a classical computer is difficult due to the exponential scaling of the Hilbert space dimension with the number of particles, it is intriguing to try and use a well-controllable quantum system with a tunable Hamiltonian to simulate the physics of another quantum system of interest. Over the last decade, this paradigm has motivated many experiments and theoretical proposals, predominantly for setups with ultracold atoms or trapped ions (see for example Ref. [4] for a recent review).
More recently, numerous proposals have emphasized the potential of the circuit QED architecture [5] to serve as an interesting solid-state quantum simulator based on polaritons [6]. Indeed, the physics and well-established fabrication techniques for coupled systems of on-chip microwave resonators and superconducting qubits, offer a number of promising properties. These include an immense freedom in engineering lattices with different geometries, dimensionalities, and topologies in a bottom-up approach [7]. Manipulation and measurement of quantum states at the single-site level are already standard practice in circuit QED experiments. Finally, the controlled coupling to local or global microwave drives also promises new insight into the nonequilibrium physics of strongly correlated polariton systems.
According to theoretical studies [8,9,10,11,12], these systems provide access to the superfluid-Mott insulator phase transition in the setting of the Jaynes-Cummings model, with characteristic similarities and differences compared to the more familiar transition in the popular Bose-Hubbard model [13]. Theoretical predictions for the physics of finite Jaynes-Cummings arrays in the presence of drive and decay exist [14,15,16,17,18,19,20,21], and proposals for introducing synthetic gauge fields [22,23] may turn Jaynes-Cummings lattices into fruitful ground for studying novel quantum phases of strongly-correlated bosons with broken time-reversal symmetry. Experimental efforts to extend the demonstration of singlesite photon blockade [24] to larger arrays are currently underway.
In this article, we study the minimal circuit QED system in which breaking of timereversal symmetry is observable: a three-site Jaynes-Cummings ring subject to a synthetic gauge field, which induces nontrivial phase factors in the photon hopping terms. Our study is organized as follows. In Section 2, we review the many-body Hamiltonian for a general circuit QED array, and derive magnitude and intrinsic signs of the photon hopping amplitudes. In Section 3, we develop the formalism to calculate linear-response reflection and transmission amplitudes. These quantities play a central role in all circuit QED experiments to date, which routinely probe quantum states by homodyne measurements. It is to be expected that transmission spectra will continue to serve as a fundamental probe in larger Jaynes-Cummings lattices, and the first experiments are currently underway ‡.
Focussing on the finite Jaynes-Cummings ring, we show in Section 4 the presence of degeneracies in the single-particle spectrum for specific values of the gauge-invariant phase sum. For commensurate polariton occupation, even weak interactions lead to strong ‡ Andrew Houck, private communication (2011). correlations and generation of highly-entangled Schrödinger cat states. In Section 5, we calculate the transmission through such a device. For weak drive strengths, we confirm transmission peaks at frequencies corresponding to the eigenmodes within the 1-excitation manifold. Beyond this linear response regime, polariton interactions become important, and numerical results show features similar to the nonlinear vacuum Rabi splitting of the singlesite Jaynes-Cummings model.

Model of circuit QED lattices
In the simplest setting possible, circuit QED arrays comprise a bosonic lattice with two species of lattice sites, cf. Fig. 1(a). The first species, qubit sites, consists of two-level systems with transition energy ε, realized by suitable superconducting qubits. The second species, resonator sites, describes harmonic oscillators with characteristic energy ω, implemented as a fixed photon mode inside on-chip resonators (e.g. coplanar waveguides).
Coupling between sites in such a circuit QED lattice allows for photon hopping and conversion between photon and qubit excitations, which in part mimics the hopping of bosonic particles in the Bose-Hubbard model [13] and its multi-species generalizations [25]. The detailed mechanisms for hopping in circuit QED lattices, however, are more diverse: hopping t between two resonator sites corresponds to the transfer of individual photons. As opposed to quantum tunneling, these processes are possible due to classical EM evanescent fields and have a clear correspondence at the classical level. For on-chip microwave resonators, this coupling is commonly realized by capacitors interconnecting two or three resonators. Hopping between resonator and qubit sites g is of the Jaynes-Cummings type, and enables the interconversion of qubit excitations and photons. Direct coupling between qubit sites should be feasible for closely spaced qubits, but, as of now, has not been demonstrated in the circuit QED architecture.
The common situation of small hopping strengths (t, g ≪ ω, ǫ) affords the rotating-wave approximation, and results in conservation of the total excitation numberN = jâ † jâ j + jσ + jσ − j , corresponding to the particle number in Bose-Hubbard type models. In total, the Hamiltonian of a circuit QED lattice assumes the form We note that disorder in the system parameters (beyond the scope of this paper), may need to be included in Eq. (1), once it is characterized experimentally. A lattice topology of particular interest, is that of a finite-size ring, shown in Fig. 1(a). Its analog for ultracold atoms, see Fig. 1(b), has been discussed in the literature [26,27,28,29], and (if loaded with bosonic atoms) shares some features with the system to be discussed here. The option to study broken time-reversal symmetry in circuit QED lattices is based on synthetic gauge fields which give rise to complex phases e iθ jj ′ in the photon hopping terms. Making photons susceptible to a magnetic field is possible, for example, by substituting interresonator coupling capacitances with passive superconducting circuits [22]. Each such circuit consists of a superconducting loop interrupted by three (or more) Josephson junctions, in which time-reversal symmetry can be broken by an external magnetic flux. Photon hopping takes place via virtual excitations of this coupler circuit, and thus transfers the symmetry breaking to the photonic system. Alternatively, active non-reciprocal devices may be used, as recently pointed out in an interesting proposal by Kamal et al [23].

Photon hopping: complex phases and signs
Individually, the phases e iθ jj ′ associated with hopping terms are not gauge invariant. Gauge transformations of the typeâ j →â j e iχ j lead to a fully equivalent Hamiltonian with altered phase factors e iθ jj ′ +iχ j ′ −iχ j . The relevant quantities are gauge-invariant phase sums, around plaquette paths C[jj ′ ] in the resonator sublattice. If no such plaquettes exist (as is the case, e.g., in an open 1d resonator chain), all phases can be gauged away immediately and the hopping terms may be rewritten with positive, real-valued amplitudes, i.e. t j,j ′ (â † jâ j ′ + H.c.). For clarity, and to make relevant signs explicit, we will assume t to be real-valued and positive throughout; i.e. phases and minus signs will not be absorbed into t.
As long as time-reversal symmetry is intact, θ Σ (C) is an integer multiple of π on all plaquettes, and hopping amplitudes can be chosen strictly real-valued. This, however, leaves open an important question about the signs of photon hopping terms. Such signs can arise from coupling λ/2 resonator modes to each other, or, in general, any half-integer mode with L = (n + 1/2)λ (n ∈ N). As illustrated in Fig. 2, it is essentially the sign change of the mode function for half-wavelength resonances which introduces a gauge-invariant phase θ Σ = π for each plaquette with an odd number of resonators. Photon hopping in resonator arrays. Photon lattices, built from quasi-1d microwave resonators, admit uniform photon hopping due to capacitive coupling between adjacent ends of two or three resonators. When coupling full-wavelength modes (resonator length L matches an integer multiple of the wavelength λ), photon hopping has a positive sign. For halfwavelength modes, L = (n + 1/2)λ, the mode function has different signs at opposite ends of each resonator, and photon hopping terms can acquire a gauge-invariant phase of π. It is important to note that the 1d nature of resonators allows two distinct configurations for coupling 3 resonators to each other, see (a) and (b). Gauge invariant π phases only emerge in ring-type coupling (b), and only for rings with an odd number of sites. Rings with even numbers of sites, see (c), do not acquire a π phase.
More specifically, photon hopping between two adjacent resonators j and j ′ has the amplitude where C c is the coupling capacitance and φ j (x) denotes the classical mode function of resonator j. For full-wavelength modes, the signs of the mode functions can be chosen to be positive at each coupling point such that hopping amplitudes are positive. Alternatively, for any different sign choice of each mode function φ j , one can always perform a gauge transformation that renders all hopping amplitudes real and positive. By contrast, half-wavelength modes are associated with an alternating sign, φ j (x end,1 ) = −φ j (x end,2 ), between opposite ends of each resonator, see Fig. 2. The additional signs lead to changes in the energy spectrum or band structure if (and only if) the lattice contains at least one plaquette with an odd number of sites on its boundary. For such odd-number plaquettes, the gauge-invariant phase sum is θ Σ (mod 2π) = π.
Derivation of hopping amplitudes.-We extend the treatment of a single, capacitively coupled resonator [22] to that appropriate for a full circuit QED lattice, and derive photon hopping amplitudes and their signs in a circuit QED lattice. Based on the usual circuit quantization scheme [30], the Lagrangian for the resonator sublattice is obtained from the LC decomposition of the transmission line resonators and a subsequent continuum limit: Here,Φ j (x, t) is the voltage profile in resonator j, and ℓ, c denote inductance and capacitance per unit length. The sum in the last term runs over nearest-neighbor resonators, C c is the coupling capacitance between each adjacent resonator pair, and the amplitudes are evaluated at the corresponding resonator ends connected to C c . Note that the sign of the coupling term stems from the expansion of the capacitive energy 1 2 C c (Φ j −Φ j ′ ) 2 | ends ; diagonal terms of the former expression appear in the Lagrangian (4), and also affect the boundary conditions for the fields Φ j [22].
Employing these boundary conditions at the two ports of every resonator, the resonator Lagrangian is diagonalized in terms of internal eigenmodes, Φ j (x, t) = µ ζ jµ (t)e iωµt φ jµ (x), and takes the form In the typical situation of weak coupling between resonators, C c ≪ C res = c L, we may focus on a single mode, and thus drop mode indices in the following.
To switch to the Hamiltonian picture, it is convenient to cast the Lagrangian into matrix | ends captures the coupling between resonators, and δ j,j ′ denotes the adjacency matrix of the resonator lattice. Legendre transform results in the Hamilton function where the vector π collects the conjugate momenta π j = ∂L/∂ζ j . The last approximation in Eq. (6) retains the leading order contribution from the resonator coupling. Finally, canonical quantization replaces ζ j → (â j +â † j )/ √ 2ω, π j → (iâ † j − iâ j ) ω/2, and in conjunction with the rotating-wave approximation, yields the resonator lattice Hamiltonian Expressed in proper frequency units, the photon hopping due to direct capacitive coupling is therefore given by where ν = ω/2π denotes the resonator frequency. As indicated in the discussion above, the signs for photon hopping are dictated by the product of mode functions involved in the hopping; specifically This result shows that, up to irrelevant sign changes induced by gauge transformations, photon hopping is strictly positive for full-wavelength modes. By contrast, for half-wavelength modes, alternating signs of the mode functions can indeed lead to observable changes. In future experiments, this may offer a convenient method for switching between photon systems with and without frustrated hopping, by a mere doubling of the microwave drive frequency.

Calculation of steady-state transmission rates
Experimentally, transmission and reflection measurements are by far the simplest probes for circuit QED arrays, and are performed routinely in small circuit QED systems with one or two resonators [31,32,33]. Conceptually, such measurements are implemented by coupling one resonator (or several) to an external transmission line, through which microwave photons can be injected, and subsequently detected in heterodyne or homodyne. For theory, however, the analysis of even this simple setting is challenging for larger numbers of resonators and qubits, since it generally involves both strong correlations, and the complexity of non-equilibrium physics for a driven, damped quantum system. Assuming that the coupling to the external transmission lines (and any additional baths) is weak, the system may be described by a Lindblad master equation of the forṁ̺ Here, the last three terms involve the damping superoperator D[L]̺ =L̺L † − {L †L ,̺}/2 [34] and account for photon loss, qubit relaxation and pure dephasing of qubits. For simplicity, we will neglect pure dephasing (γ ϕ = 0) and assume decoherence with uniform rates induced by coupling to separate baths. Extensions to shared baths and varying decoherence rates are possible, but can add further complexity to the analysis. The unitary evolution is captured by the Liouville term, and expressed in terms of the system Hamiltonian in the frame co-rotating with the drive, i.e.
Here, the unitary isÛ = e iω d tN and removes the time-dependence of the drive Hamiltonian for a single drive frequency ω d , andĤ =Ĥ lat +Ĥ d is the Hamiltonian including the drive. Under steady-state conditions,̺ = 0, the master equation (10) forms a linear system of equations for the matrix elements of̺. Its solution, in principle, enables the prediction of transmission and reflection signals by calculating â = tr(â̺), whose real and imaginary part correspond to the two quadratures of the homodyne signal in transmission. It is important to stress, however, that the brute-force numerical solution of this equation represents a formidable challenge even for arrays with only a handful of resonators and qubits: the number of independent real variables to be determined from̺ = 0 is given by (Λ + 1) 2Nr 2 2Nq − 1, where N r and N q denote the number of resonators and qubits in the lattice, and Λ is the cutoff in photon number. As an example: a lattice consisting of 3 resonators and 3 qubits leads to a linear system of equations with 10 5 − 1 variables when admitting up to 4 photons in each resonator. A more appropriate cutoff scheme, not based on a truncation of the photon Hilbert space, will be described in Section 5.
In light of this scaling, it is useful to note that the direct calculation of the full reduced steady-state density matrix can be avoided, and instead, only those elements of̺ be evaluated that indeed contribute to the expectation value â . To see this, consider the vanishing time derivative of the steady-state expectation value, and re-express it via the master equation: Not surprisingly, this equation for a j in general does not close after explicit evaluation of the commutator. For interacting systems it instead produces an infinite hierarchy of equations, involving an ever increasing number of expectation values of additional operators, which is analogous to situation of the equations of motion method for the system's Green's functions. We review the simple and exact solution in linear response first, and discuss numerical results beyond this in Section 5. Interaction in circuit QED lattices is readily understood as a hard-core repulsion between bosonic excitations on qubit sites, which is so strong that levels beyond the ground and first excited state are pushed to inaccessibly high energies. Formally, the Jaynes-Cummings lattice thus corresponds to a two-species Bose-Hubbard model with U 0 = 0 on resonator sites and U 1 → ∞ on 2-level sites [35]. This interaction, of course, is only effective in the presence of multiple "particles" (i.e. polariton excitations). Hence, it may be ignored in the linear response regime where weak driving does not lead to multiple excitations. In this limit, all qubit sites can be replaced by a second set of harmonic oscillator sites,σ − j →b j , and the resulting lattice Hamiltonian, describes a simple system of coupled oscillators. Using ∆ υ = υ −ω d to denote the detuning of the drive from frequency υ, the commutator required for the calculation of â j , Eq. (13), produces where nn(j) is the set of sites that are nearest neighbor to site j. Together with the analogous relations for d b j /dt, one obtains the matrix equation Here, we used the following notation: homodyne amplitudes for the two kinds of sites are lumped into the vectors α = ( â 1 , . . . , â n ) and β = ( b 1 , . . . , b n ′ ). Likewise, Ω = (Ω 1 , . . . , Ω n ) collects the drive strength for each resonator (â j sites). Finally, T is the resonator coupling matrix defined by T jj ′ = te iθ jj ′ if j ′ ∈ nn(j), and = 0 otherwise.
From Eq. (16), one obtains the exact solution for the resonator homodyne signals as For the physical interpretation of this expression, it is appropriate to switch to the eigenmode basis of the resonator subsystem, U(∆ ω + T)U † = diag(∆ Υν ), where ∆ Υν denotes the drive detuning from eigenmode frequency Υ ν , and U † is the matrix of eigenmode amplitudes for the resonators; specifically, U * νj = ( ε ν ) j gives the relative amplitude of eigenmode ν in resonator j. Using this basis yields the following expression for the complex homodyne amplitudes The general structure of this expression can be explained as follows. Coherent driving of individual resonators (with amplitudes given by the components of Ω) results in an effective drive Ω Υ = U Ω of the individual resonator eigenmodes. The complex amplitudes induced in each eigenmode ν depend on the detuning of the drive frequency from photon modes and qubits as well as the dissipative rates for photon loss and qubit relaxation. Finally, the amplitude measured in resonator j is the weighted sum over all eigenmode amplitudes, with weight factors U * νj accounting for the overlap between resonator j and eigenmode ν.

Resonator array (g = 0)
For g = 0, transmission only probes the array of coupled resonators. As expected, a drive on a single resonator k results in a transmitted amplitude T jk (ω d ) = | â j | in resonator j, which becomes maximal whenever ω d matches one of the eigenfrequencies Υ ν . The maximum amplitude values (peak heights in the homodyne amplitude) are given by where Ω (ν) jk = |U * νj U νk Ω k | denotes the effective drive strength on resonator j, mediated by eigenmode ν, and generated by a drive on resonator k. The final approximation of Eq. (19), is valid for non-degenerate and well-separated eigenmodes, i.e. |Υ ν ′ − Υ ν | ∼ t ≫ κ. The generalization to degenerate but otherwise well-separated eigenspaces is achieved by including an additional sum over those modes degenerate with ν.

Jaynes-Cummings lattice (g = 0)
Once the qubit subsystem is coupled to the resonator array, the single-excitation spectrum and corresponding transmission resonances are modified. Focussing again on sufficient peak separation (g, t ≫ κ, γ), the peak frequencies can be approximated by the poles of Eq. (18) when setting κ = γ = 0. This yields as the predicted resonance frequencies for the transmitted amplitude. It is simple to verify that the frequencies Θ ν,± are identical to the eigenfrequencies of the full coupled system of the resonatorsâ j and the harmonic modesb j . The corresponding heights of resonances in the homodyne transmission amplitude are The Jaynes-Cummings physics of resonant vacuum Rabi splitting, and off-resonant dispersive shifts evidently carries over to the single-excitation physics of the Jaynes-Cummings lattice. The one-excitation spectrum, described by the set of frequencies {Θ ν,± }, crucially depends on the detunings of the qubit frequency ǫ from each resonator eigenmode Υ ν . For each individual eigenmode ν, detuning can be large, |ǫ − Υ ν | ≫ g, and render the interaction dispersive, or be near resonant, |ǫ − Υ ν | ≪ g, and lead to strong photon-qubit hybridization.
Dispersive regime.-The dispersive coupling between qubits and resonator mode ν at large detuning is reflected in small frequency shifts of the qubit frequency ǫ and the eigenmode Υ ν . Expansion of Eq. (20) to leading order in g/(ǫ − Υ ν ) yields the dispersive shifts familiar from the Jaynes-Cummings model, which imply level repulsion between photon and qubit modes. The corresponding height of the photon resonance is only slightly reduced from its g = 0 value, whereas the additional qubit resonance is strongly suppressed: Resonant regime.-Whenever the qubit frequency is tuned into near degeneracy with one of the photon eigenmodes, their coupling results in hybridized eigenstates and the typical vacuum Rabi splitting. Expansion of equations (20) and (21) to leading order in (ǫ − Υ ν )/g gives for the peak frequencies and heights. It is interesting to note that this linearization approach, valid for circuit QED arrays at sufficiently small drive powers, may also provide insight beyond the linear response regime for multi-level systems with small anharmonicity. The transmon qubit [32,36] provides one example for such a system, where under appropriate conditions the anharmonicity may be treated perturbatively [15].

Strongly-correlated states in a three-site Jaynes-Cummings ring lattice
The discussion in Sections 2 and 3 has been general. In the following, we will focus on the concrete example of a three-site Jaynes-Cummings ring lattice, see Fig. 3(a). The system consists of three resonators coupled via a Josephson ring circuit [22], and each resonator interacts with its local superconducting qubit. It is the simplest system for which a gauge field can lead to measurable effects as the nonzero phase sum around the loop cannot be gauged Josephson ring circuit, and each resonator interacts with its local superconducting qubit. In Section 5 we will consider the case where one of the resonators is driven by a microwave tone. (b) Spectrum in the 1-excitation manifold for g = 0 (solid) and g/t = 0.2 (dashed). For θ = π/3 the highest-energy state is doubly degenerate. (c) Spectrum in the 3-excitation manifold. The degeneracy in the 1-excitation manifold at θ = π/3 induces a 4-fold degeneracy in the 3-excitation manifold which is lifted for g/t = 0.5. (d) Overlap of the highest-energy 3-excitation state with |3, 0, 0 | ↓, ↓, ↓ (red solid) and |0, 3, 0 | ↓, ↓, ↓ (green dashed) as a function of phase twist θ for g/t = 0.5. Around θ = π/3 the state has large overlap with the superposition states |ψ + in Eq. (29). The insets in (c) and (d) show zoom-outs of the same graphs. In all cases, the resonators and qubits are in resonance with each other ε = ω.
away. In this section we show how strongly-correlated states can arise even in this simple setup.
Following Section 2, the Hamiltonian of the three-site Jaynes-Cummings ring lattice is given bŷ where θ = θ Σ /3 denotes the phase twist and θ Σ is the gauge-invariant phase sum. We note that a similar Hamiltonian arises for bosons in rotating ring lattices [27,28,29]. SinceĤ 0 (25) conserves the total excitation numberN = 3 j=1 â † jâ j +σ + jσ − j , its energy spectrum and eigenfunctions can be studied for each N separately. In the subspace with one excitation (N = 1) we can regard the qubits as a second class of harmonic oscillators: σ j →b j . Since the system has translational symmetry it is advantageous to introduce quasimomentum operatorsâ j = 1 √ 3 2 k=0 e −2πikj/3â k andb j = 1 √ 3 2 k=0 e −2πikj/3b k . These reduce the Hamiltonian to the form with the eigenmode frequencies Υ k = ω+2t cos 2πk 3 + θ of the coupled resonator sublattice, and eigenvalues Fig. 3(b) shows the 1-excitation energy spectrum as a function of the phase twist θ. For g = 0, the spectrum consists of the resonator modes Υ k and the three uncoupled qubits at ǫ = ω. For g = 0 several level crossings between qubit and photon excitations turn into avoided crossings, but at the critical phase twist θ = 0, the ground state (and for θ = π/3, the highest-energy state) in the 1-excitation manifold remains doubly degenerate, i.e. Θ − 1 = Θ − 2 for θ = 0 and Θ + 0 = Θ + 1 for θ = π/3. Accordingly, even without synthetic gauge field, θ = 0, the plaquette is frustrated. The degeneracy in the 1-excitation spectrum in turn leads to strong correlations among the photons even for weak interaction strength, as we will show next.
If (and only if) the number of lattice sites L and the number of excitations N are commensurate, these degeneracies are lifted by the Jaynes-Cummings term for g = 0. To understand the induced correlations, we rewrite the Hamiltonian (25) aŝ The last term lifts the degeneracy between |N, 0, 0 | ↓, ↓, ↓ and |0, N, 0 | ↓, ↓, ↓ in a 2Norder process, and the new highest-energy eigenstates are the cat-like superposition states
We see that synthetic gauge fields can lead to degeneracies in the 1-excitation spectrum which may induce strong correlations in the interacting many-excitation system. | 2 , in response to driving at port 1, | â in 1 | 2 : in (a) for g = 0 as a function of drive detuning for ∆ ω for θ = 0 (solid) and θ = π/6 (dashed) and in (d) for g/κ = 1 as well as δ/t = 1 (solid) and δ/t = 2 (dashed). Transmission for g = 0 as a function of phase twist θ and drive detuning ∆ ω is shown in (b) and for g/κ = 1 as a function of detuning δ and drive detuning ∆ ω in (e). For comparison, we plot the spectrum of the 1-excitation manifold as a function of phase twist θ for δ = g = 0 in (c) and as a function of detuning δ for g/t = 1/3 and θ = 0 in (f). The other parameters are t/κ = 3 and γ/κ = 0.1. In (c) and (f) resonator-like and qubit-like eigenstates are indicated by solid and dashed lines, respectively.

Transmission of a three-site Jaynes-Cummings ring lattice
To establish contact with experimental measurements of homodyne transmission, we apply the general discussion of Section 3 to the finite Jaynes-Cummings ring. As indicated in Fig. 3(a), we consider driving a single resonator with a microwave tone of strength Ω. All resonators and qubits are subject to dissipation leading to line widths κ and γ, respectively.
As mentioned above, the eigenmodes of the resonator array are characterized by the quasi-momentum number k, and the amplitude of eigenmode k in resonator j is U * kj = ( ε k ) j = e −2πikj/L / √ L, independent of phase twist θ. Using the input-output relationŝ a out j =â in j + √ κ jâj and √ κ j â in j = iΩ j [37], we obtain the reflection and transmission rates from the expectation values of the photon annihilation operators, â j . In reflection, the drive interferes with the reflected wave from the cavity, i.e.
In the absence of coupling between resonator and qubits (g = 0) or if qubit dissipation is negligible (γ = 0), we have j n out j = n in 1 expressing the conservation of the photon flux n in/out j = (â † j ) in/out (â j ) in/out . We first consider the case g = 0, where qubits and resonators are decoupled. Fig. 4(a) shows the transmission | â 2 | 2 as a function of detuning between resonator and drive frequency ∆ ω . As expected, peaks in transmission correspond to the eigenfrequencies of the resonator array. For θ = 0 the quasi-momentum states k = ±1 are degenerate leading to only two peaks at −∆ ω /t = −1 and −∆ ω /t = +2. For θ = 0 the spectrum is in general nondegenerate and the three peaks mark the three quasi-momentum modes of the ring lattice. In Fig. 4(b), the dependence of transmission on drive detuning ∆ ω and phase twist θ shows that the transmission follows the resonator dispersion ε k (θ), cf. Fig. 4(c). An interesting detail: when the spectrum is almost degenerate transmission is very sensitive to changes in the phase twist, see, e.g., the situation of −∆ ω /t = 1 close to θ = π/3. For these parameters the circuit, in fact, implements an on-chip circulator (different from the parameter choice of Ref. [22]).
For g = 0 and small drive, linear response theory remains valid and qubit sites can be replaced by harmonic oscillators as discussed in Section 3. Results from this treatment are shown in Fig. 4(d) where the transmission | â 2 | 2 is plotted as a function of detuning ∆ ω (g/t = 1 and fixed θ = 0). For δ = ǫ − ω = t, the qubits are not in resonance with any of the eigenmodes of the resonator array. Relative to the non-interacting case g = 0, the coupling to the qubits slightly shifts the frequency of the resonator peaks and three additional peaks (two of them degenerate at θ = 0) appear in the spectrum close to −∆ ω = δ. By contrast, for δ/t = 2 the qubits are in resonance with the k = 0 resonator mode and, analogous to the vacuum Rabi splitting in the single-site Jaynes-Cummings model, the resonance peak splits. In Fig. 4(e) we show the transmission | â 2 | 2 as a function of the detuning δ between qubit and resonator and drive detuning ∆ ω . As the qubit frequency ǫ is varied, qubits are tuned in and out of resonance with each of the resonator modes. For comparison, Fig. 4(f) depicts the spectrum of the system in the 1-excitation manifold as a function of qubit detuning δ.
For stronger drive Ω, population beyond the 1-excitation subspace sets in, and linear response theory breaks down. To find the transmission in that case, we compute the relevant expectation values via the steady-state solution of the full quantum master equation (10). As mentioned above, even small resonator arrays and moderate photon number cutoffs lead to matrix sizes for the superoperator which quickly become impractical. This problem can be mitigated slightly by an adaptive truncation scheme for the density matrix elements. First, note that the HamiltonianĤ 0 conserves the number of excitations N, while the drive HamiltonianĤ d induces transitions between subspaces of different excitation numbers. Each element of the standard Fock basis |µ is associated with an excitation number N. Now, assigning to each element of the density matrix ̺ µν = µ|̺|ν a generalized excitation number N defined as the sum of the excitation numbers of the states |µ and |ν , a truncation based onÑ is possible where the cutoff is adjusted with increasing drive strength Ω. This cutoff scheme is more appropriate and fine-grained compared to a maximum photon number cutoff in Fock space.
In Fig. 5(a) and (b) we present the transmission | â 2 | 2 as a function of drive detuning ∆ ω for different drive strengths Ω/κ (cutoff:Ñ = 4). For weak drive, the response is close to the , Ω/κ = 0.5 (green dash-dotted), and Ω/κ = 1 (red dashed). Other parameters are t/κ = 5, g/κ = 100, δ/κ = 5, and γ = θ = 0. Blue and red vertical lines indicate the one-photon and one half of the two-photon transition energies, respectively. (c) For various driving strengths, probability P 0 (blue), P 1 (green), and P 2 (red) to have zero, one, and two excitations in the system. Inset in (a) illustrates the fact that the 2-photon transition energy |0 → |2 is detuned from twice the 1-photon transition |0 → |1 . linear-response result. As the drive increases, the transmission is suppressed compared to the linear response value. As better seen in Fig. 5(b), the peak additionally becomes asymmetric, shifts to smaller detunings and finally at strong drive splits into two. In Fig. 5 we have marked the 1-photon transition energies and half of the 2-photon transition energies, revealing that the second peak appears at half the 2-photon transition frequency. This is further corroborated in Fig. 5(c) where the probabilities for zero, one and two excitations are plotted. At the frequency of the 1-photon transition increasingly strong driving leads to a larger population of the 1-excitation sector P 1 but the probability for two excitations P 2 remains small. At half the 2-photon transition frequency the population of the 1-excitation space P 1 remains small, but at stronger drive strength the 2-excitation manifold is occupied.
This behavior resembles the nonlinear vacuum Rabi splitting as observed in a singlesite Jaynes-Cummings system [38]. Since the spectrum is nonlinear, photon blockade [24] prevents population transfer from the 1-excitation to the 2-excitation subspace due to the energy mismatch of E 1 − E 0 vs. E 2 − E 1 , see inset of Fig.5(a). Specifically, driving the system on the 1-photon transition |0 → |1 transfers population from the ground state |0 to the first excited state |1 , but a second excitation with the same energy cannot enter the system. On the other hand, driving the system at half the 2-photon transition frequency |0 → |2 , not very effective at low power, allows for two-photon transitions at sufficiently large power. The features of the transmission signal can thus be explained naturally as a consequence of the nonlinear spectrum of the system.

Conclusion
In summary, we have explored the physics of a Jaynes-Cummings ring, the minimal circuit QED system manifesting effects of broken time-reversal symmetry. Analogous to bosons in rotating ring lattices, degeneracies at specific values of the gauge field and interactions lead to strongly-correlated states. As future experiments are likely to employ transmission as one probe for such systems, we have calculated the response to drive and decay, generalizing phenomena known from the single-site Jaynes-Cummings model, e.g. the vacuum Rabi and nonlinear vacuum Rabi splittings. It is clear that a theoretical description of transmission experiments with larger lattices and/or higher drive power will require novel simulation techniques beyond the explicit calculation of the steady-state density matrix of the master equation. Analytical and numerical methods from the field of quantum transport may turn out to be adaptable to this novel setting. Other valuable experimental probes beyond transmission, such as two-tone spectroscopy and correlation functions, reveal additional information about the system, and will be explored in future work.