Bose-Hubbard dynamics of polaritons in a chain of circuit QED cavities

We investigate a chain of superconducting stripline resonators, each interacting with a transmon qubit, that are capacitively coupled in a row. We show that the dynamics of this system can be described by a Bose-Hubbard Hamiltonian with attractive interactions for polaritons, superpositions of photons and qubit excitations. This setup we envisage constitutes one of the first platforms where all technological components that are needed to experimentally study chains of strongly interacting polaritons have already been realized. By driving the first stripline resonator with a microwave source and detecting the output field of the last stripline resonator one can spectroscopically probe properties of the system in the driven dissipative regime. We calculate the stationary polariton density and density-density correlations $g^{(2)}$ for the last cavity which can be measured via the output field. Our results display a transition from a coherent to a quantum field as the ratio of on site interactions to driving strength is increased.


Introduction
In recent years, the investigation of condensed matter and quantum many-body systems with quantum simulators, artificial quantum many-body systems that offer unprecedented controllability and measurement access in the laboratory, has become an active research area and is currently receiving increasing attention. The technology employed for quantum simulators ranges from ultra-cold atoms [1] to ion traps [2] and systems of coupled cavities [3,4,5,6] among others.
Due to recent technological progress, arrays of coupled cavities and optical nano-fibers in which the trapped light modes couple to atoms are now becoming suitable devices for the generation of quantum many-body systems of polaritons [3,4,5,6,7,8,9,10], i.e. quantum-mechanical superpositions of atomic and photonic excitations. In these systems, it is of particular interest but also most challenging to reach a strongly correlated regime, where their dynamics differs most significantly from that of classical light fields. The key experimental requirement for reaching these conditions is a so called strong coupling regime for the cavities respectively the fiber. This means that the coherent coupling between the light modes and the atoms or other optical emitters must be strong compared to the loss processes which are inevitably present in every device. A very impressive strong coupling regime has recently been realized in circuit cavities, making these devices an ideal platform for studying strongly correlated polaritons.
Circuit QED [13,14,15] was developed as a solid state equivalent to optical cavity QED, coupling Josephson qubits that are acting as artificial atoms with stripline resonators acting as cavities for microwave photons. Here, the reduced quasi onedimensional mode volume of the stripline resonator and the enhanced dipole moment of the Josephson qubit with respect to atoms give rise to a pronounced strong coupling regime, where the coupling between resonator and qubit, g, significantly exceeds both the decay rate of the resonator, κ, and the qubit, γ, g/κ 1 and g/γ 1. Moreover, since stripline resonators trap microwave photons, they are more than 1cm long. The precision of current fabrication techniques thus allows to build several resonators that can resonantly couple to each other via mutual photon tunneling on the same chip [16,11,12]. The currently employed circuit QED technology thus permits to build arrays of resonantly coupled cavities that each interact in a strong coupling regime with qubits. In this way it is one of the first setups to feature all properties which are needed to generate strongly correlated many-body systems of polaritons.
In this work we show that an effective Bose-Hubbard Hamiltonian for polaritons can be engineered in an array of stripline resonators that each couple to a transmon qubit [17]. Josephson qubits [18,19] come in basically three different flavours depending on the property that is controlled from the outside or rather the channel of the qubit environment coupling: flux- [20], phase- [21] and charge-qubits [22,23]. Here we consider a setup with transmon qubits [17] which are charge qubits (Cooper pair boxes) operated at sufficiently enhanced values for the ratio of Josephson energy, E J , over charging energy, E C , E J /E C ≥ 50 and are robust against decoherence caused by fluctuations of background charges. We emphasize that the effective Bose-Hubbard Hamiltonian we derive can be realized with resonators and qubits of readily existing technology.
Quantum phases for the ground state and low temperature thermal states of the Bose-Hubbard Hamiltonian have been studied with ultra-cold atoms trapped in optical lattices [1]. This system has also been employed to study the dynamics of noneequilibrium states that were prepared by sudden quenches of some lattice parameters [24]. In contrast, a realisation in an array of stripline resonators allows to investigate the Bose-Hubbard Hamiltonian in a fundamentally different regime, where the resonator array is permanently driven by lasers to load it with photons and thus compensate for the excitations that are lost due to qubit relaxation and cavity decay. Whereas substantial understanding of equilibrium quantum phase transitions has been achieved, a lot less is known about these non-equilibrium scenarios where the dynamical balance between loading and loss mechanisms leads to stationary states. It is the investigation of these stationary states, that our approach to the Bose-Hubbard Hamiltonian is ideally suited for.
Experiments with transmon qubits [17,25] coupled to a stripline resonator are often conducted without directly measuring the state of the qubit but by spectroscopically probing the transmission properties of the resonator. In an experiment the effective Bose-Hubbard Hamiltonian will thus be operated out of thermal equilibrium in a driven dissipative regime [26,27,28,29]. In a suitable setup with a linear chain of resonators one would thus drive the first resonator with a coherent microwave input and measure the properties of the output signal at the opposite end of the chain. Figure 1. Sketch of the proposed system to simulate Bose-Hubbard physics. Stripline resonators are coupled capacitively in a chain and each stripline resonator is coupled to a transmon qubit that gives rise to an on-site interaction for the polaritons. For definitions of J 0 and g, see (7) and (4).
In the regime we consider, this situation can be accurately described by a Bose-Hubbard Hamiltonian for polaritons with a coherent driving term at the first site and Markovian losses of polaritons due to cavity decay and qubit relaxation at all sites of the chain. In this scenario, the interplay of coherent drive and polariton loss leads to the emergence of steady states, for which we derive the particle statistics and characteristic correlations. In doing so we focus on the polariton statistics, in particular the density and density-density correlations, in the last resonator as these can be measured via the output signal. Our results show a transition from a coherent field to a field with strongly non-classical particle statistics as the ratio of on site interactions to driving strength is increased.
This work is devided into two main parts. In section 2 we show that the dynamics of our system can be described by a Bose-Hubbard model and in section 3 we present the results of our calculations for the polariton density and density-density correlations.

Transmon-QED and the Bose-Hubbard Model
To generate a Bose-Hubbard model with polaritons we consider an array of capacitively coupled stripline resonators with each resonator coupled to a transmon qubit, see figure 1. In this section, we first introduce the Hamiltonian that describes this setup and then show how it can be considerably simplified and transformed into a Bose-Hubbard Hamiltonian for two polariton species.

The full Hamiltonian
The full Hamiltonian of our setup is a sum of single-site Hamiltonians, H 1−site,i , that each describe a transmon qubit coupled to a stripline resonator and terms that describe the capacitive coupling between neighbouring stripline resonators, H C J ,i,i+1 , A transmon qubit can be regarded as a Cooper pair box that is operated at a large ratio of E J /E C 1. This regime can be accessed by shunting the Josephson junction with an additional large capacitance and thereby lowering the charging energy E C = e 2 /(2C Σ ).
Here, C Σ = C J + C g + C B is the sum of the junction's capacitance, C J , the mutual capacitance with stripline resonator, C g , and the shunting capacitance C B . The Hamiltonian for one stripline resonator coupled to a transmon qubit reads, Here, ω r is the resonance frequency of the isolated resonator and we have omitted the site-index i for readability. Transmon qubits consist of two superconducting islands connected by Josephson junctions.n is the operator for the difference in the number of Cooper pairs on the two superconducting islands, n dc g the offset charge induced by an applied dc voltage and intrinsical defects andφ is the operator for the superconductingphase difference between the two islands. We assume the transmon qubit to be placed in the antinode of the stripline resonator's field mode. This gives rise to an additional ac component in the offset charge, n ac rms = ω r /C r is the root mean square voltage of the vacuum field mode and a the annihilation operator of photons in the stripline resonator. The offset charge n ac g thus induces a coupling between the transmon and photons in the resonator. For circuit QED setups one normally uses λ-resonators with the antinode located at the middle of the resonator.
The energy of the coupling capacitor between neighbouring transmission line resonators, e.g. sites i and i + 1, can be expressed in terms of the difference in the electrostatic potentials across the capacitor, Here, C r is the capacitance of the whole stripline resonator with respect to the ground plane and C J the capacitance of the capacitor that connects the two resonators. We assume the electrostatic potential in resonator i to have antinodes at the ends of the resonator and write it in terms of the creation and annihilation operators, a † i and a i . We now turn to simplify the Hamiltonian (1) by a sequence of approximations.

Approximations to single-site terms
We first simplify the single-site terms, H 1−site , as in (2). For large E J /E C and low energies, the phase difference between the two islands remains small and we can expand the cosine in (2) around ϕ = 0 up to quartic order, Higher order terms can be neglected, c.f. [17]. In terms of bosonic creation and annihilation operators for the transmon qubit excitations, where with, The terms linear in the creation and annihilation operators can be eliminated by performing the unitary transformation that displaces the creation and annihilation operators by the constants r and s respectively is. r and s can now be chosen such that all terms linear in a and b cancel in the transformed Hamiltonian. Finally the interaction between the transmon qubit and the field mode of the stripline resonator gets reduced to an exchange interaction in a rotating wave approximation. To justify this rotating wave approximation we have to ensure that the interaction strength between the transmon qubit and the stripline resonator is smaller than the sum of the frequencies of the two, g ω r + ω q 1 .
Parameters extracted from [30] are ω r = 43.6Ghz, E C = 0.4Ghz and a maximal value for E J /E C of 150. We choose C g eV rms /(C Σ ω r ) = 0.1 which is in agreement with the theoretical upper bound in [17] and find g ω r +ωq ≈ 0.1. The single-site Hamiltonian can thus be approximated by,

Approximations to couplings between resonators
We now turn to simplify the couplings between neighbouring resonators, H C J ,i,i+1 . We assume that C J C r which implies that C J ω r /(2C r ) is small compared to the isolated cavity frequency ω r , C J /(2C r ) 1 and apply a rotating wave approximation to neglect those terms in the intercavity interaction that don't conserve the total photon number, where J 0 = C J Cr ω r . The first term on the rhs of (7) can be absorbed into the single-site Hamiltonians by introducing a shifted resonator frequency and the remaining term in (7) describes tunneling of photons between neighbouring resonators. Next, we explain how the simplified Hamiltonian can be transformed to a two component Bose-Hubbard Hamiltonian.

The polariton modes
In the case of circuit QED with transmon qubits the coupling constant between photonic and qubit excitations is the dominating interaction energy of the system. Excitations of the whole system therefore can't be characterized as purely photonic or qubit excitations in general. To obtain a more suitable description we introduce new creation and annihilation operators, describing excitations commonly termed polaritons where, sin(θ) = g g 2 + ∆ω + ∆ω 2 + g 2 2 cos(θ) = ∆ω + ∆ω 2 + g 2 with ∆ω = ω r − ω q . The sine and cosine terms account for the transition of the character of the exitations from photonic to qubit excitations for the c + -mode as the ratio E J /E C increases and vice versa for the c − -mode. Expressing the Hamiltonian C J ,i,i+1 in the polariton modes (9a-b) we get, This Hamiltonian consists of two harmonic chains for the c + and c − polariton modes, with ω ± = ((ω r + ω q ) ± ∆ 2 + g 2 )/2, a term describing hopping from a c − -mode at site i to a c + -mode at site i + 1 and all other possible combinations, and a term describing the nonlinearity, We assume the frequencies of the two polariton modes to be well separated, apply another rotating wave approximation where we neglect the term H cc and convert the nonlinearity term H nlin into Kerr form and get a renormalization of the polariton frequency and a density-density coupling between the polariton modes. This requires the difference in frequencies for the unperturbed modes ω + − ω − = ∆ω 2 + g 2 involved to exceed the magnitude of the coupling between the modes and the nonlinearity, Plugging in realistic values for the parameters, extracted for example from [30], (E C = 0.4Ghz, ω r = 43, 6Ghz E J /E C = 0..150) we realize that the second inequality is indeed fulfilled. Engineering the capacitance C J such that J 0 is of the order of E C /12, we can ensure that the first equality is fulfilled as well. The rotating wave approximation eliminates the intermode exchange coupling and we obtain a Bose-Hubbard Hamiltonian for both modes, c + and c − , with a density-density coupling between them, where with We thus arrived at a two component Bose-Hubbard model for the modes c + and c − with attractive interactions and a density-density coupling between both species. The two species are a mixture of stripline resonator field mode and qubit excitations (9a-b) with different weights of the photonic or qubit contribution depending on the value of For small values of E J /E C the c + polaritons become increasingly photonic. Consequently, their tunneling rate J + approaches the tunneling rate of bare photons, J 0 , and their on-site interaction U + vanishes. For large E J /E C , on the other hand, J + vanishes and the nonlinearity U + approaches the nonlinearity of the qubits, E C . For the c − polaritons, the roles of both limits are interchanged.
For each value of E J /E C , the separation between the resonance frequencies of c + and c − polaritons, |ω + − ω − |, is sufficiently large such that, in a scenario where we drive the first cavity by a microwave source, we can always adjust the frequency of the drive to only selectively excite one of the modes. For reasons that will become clear later we therefore choose the c + -polaritons to be our quantum simulator for a driven dissipative Bose-Hubbard model.

Validity of the approximations
To further illustrate the validity of our approximations we compare the eigenenergies of the full Hamiltonian H, c.f. (1), approximated under assumption (5), with the eigenenergies of the Bose-Hubbard Hamiltonian H (3) , c.f. (10). The single-site Hamiltonians summed up in the full Hamiltonian describe the interaction between transmon qubit and stripline resonator in a rotating wave approximation. This Hamiltonian has already been used to describe an experiment revealing the nonlinear response of a resonator and transmon qubit system with excellent agreement between theory and experimental data [33]. Therefore comparison of the eigenvalues of our Bose-Hubbard Hamiltonian and the eigenvalues of the full Hamiltonian provides a good means to estimate the effects of the approximations we made. For simplicity we restricted our model to two sites.
Both Hamiltonians conserve the total number of excitations and we can diagonalize them in each subspace with a fixed number of excitations independently. Eigenvalues of the full Hamiltonian in the one excitation subspace are plotted in solid lines in figure 3a. Without qubits, the Hamiltonian of the two resonators has eigenmodes a ± = (a 1 ± a 2 )/ √ 2. In figure 3 a) we also plotted the energies of these eigenmodes of the two coupled empty stripline resonators marked by two horizontal dash-dotted gray lines and the eigenenergy of the transmon qubit marked by a dash-dotted gray line. Additionally to the differences between the eigenenergies of the full Hamiltonian and the Bose-Hubbard Hamiltonian in the one excitation subspace we plotted the differences in the two excitation subspace in figure 3 c). These eigenenergies can be grouped for the Bose-Hubbard Hamiltonian according to the distribution of excitations among the two polariton species. Differences of eigenenergies for states with two c − polaritons are plotted in blue, for two c + polaritons in red and for one c − polariton and one c + polariton in green. In the two excitation subspace we have similar findings as in the single excitation subspace. There are aberrations for the anti-crossing area because of the neglected intermode polariton exchange interaction. In addition, states containing c − polaritons have aberrations for small values of E J /E C due to the approximations of the transmon Hamiltonian whereas c + polaritons do not.
Therefore the Bose-Hubbard Hamiltonian for the c + polaritons mimics the behaviour of the full Hamiltonian for the full range of E J /E C , provided the intersite coupling J 0 is at most of the order of the on-site nonlinearity E C and the polariton densities are not to high. To conclude: In a driven dissipative setup where we selectively excite the c + -Polaritons we do have a quantum simulator for a Bose-Hubbard Hamiltonian.

Polariton statistics in the driven dissipative regime
In this section we make use of the above explained mapping of the full Hamiltonian H to a two component Bose-Hubbard Hamiltonian H (3) and consider a chain of coupled resonators, where we coherently drive the first resonator and adjust the microwave drive frequency to selectively excite the c + -polaritons. In the driven dissipative regime we expect to explore new physics that go beyond the equilibrium features that are commonly examined in many body physics. We thus calculate the polariton density and the density-density correlations g (2) in a master equation approach and analyse the dependencies on the system parameters J + ,U + and the Rabi frequency of the microwave drive Ω.
First experimental realisations of coupled stripline resonators are expected to consist of only a few resonators. To closely approximate the expected experiments and to speed up numerical calculations, we thus focus on a minimal chain of only two resonators. More specifically, we consider two stripline resonators coupled to transmon qubits, where the first stripline resonator is driven by a microwave source and the output signal of the second cavity is monitored as a function of the microwave drive frequency and the ratio of E J /E C which can be controlled by applying an external magnetic flux to the transmon qubits c.f. figure 4. This setup and very similar setups are currently investigated in experiments for example [16], and the spectroscopic measurement technique proposed here has already been demonstrated in single site experiments for example in [14].
The ouput fields are linear functions of the intra-cavity field in the second resonator and thus show the same particle statistics. We therefore calculated the polariton density and the g (2) -function for the second cavity. To do this, we use a master equation   approach in which each element, the stripline resonators and the transmon qubits, couple to separate environments with decay rates denoted κ for the stripline resonators and γ for the transmon qubits. Absolute values can for example be extracted from [33] where γ = 3.7MHz . Decay of the stripline resonator is due to the finite transparency of the coupling capacitors at both ends of the resonators and decay rates for example in [30] are κ = 5.7MHz. Both environments, for the transmon qubit and the stripline resonator, are assumed to be in a vaccum state which is a valid assumption at typical temperatures for circuit QED experiments of T = 15mK. Therefore in a master equation for a Hamiltonian expressed in the operators for the resonator field mode a and the transmon qubit b the dissipators read, These can be cast into dissipators expressed in the polariton modes c + and c − , In the driven dissipative case where we selectively excite the polariton c + mode we can We truncate this set of coupled equations by omitting couplings to mean values with n + m + k + l bigger than some n max and solve the reduced set of equations of motion.
To confirm the accuracy of our approach, we test its convergence with increasing n max . That is, we repeat the procedure for n max → n max + 1, compare the results and increase  Figure 5. Logarithmic density plot of the polariton density c † 2 c 2 in the last cavity plotted against E J /E C and the frequency of the microwave drive in units of the frequency of the stripline resonator ω µw /ω r . Resonances in the density of polaritons arise where the microwave frequency matches one of the transition frequencies of the non-driven conservative system Hamiltonian H c+ the value for n max in case both results differ by more than some required threshold value. The advantage with respect to a method that truncates the Hilbert space at some maximal number of excitations, is that our method becomes exact in the limit where the Hamiltonian becomes harmonic which is the case for small values of E J /E C . Moreover we experience a substantial decrease in cpu-time for this method.

Polariton density
We are interested in the field particle statistics in the driven dissipative regime and its dependencies of the on-site nonlinearity U + , the intersite coupling J + and the strength of the microwave drive Ω. We therefore first consider the density of c − polaritons in the last resonator. Figure 5 shows the density of c + polaritons, c † 2 c 2 , in the second cavity as a function of the ratio E J /E C and the microwave drive frequency ω µw . The density of polaritons in the last cavity exhibits resonances when the microwave drive frequency matches one of the transition energies of the undriven conservative system Hamiltonian H c + and decreases rapidly because of the small decay rate Γ. One can clearly see the resonances due to transitions driven between the groundstate and eigenenergies in the one excitation subspace plotted in figure 3 a). Transitions from the groundstate into a two excitation state are much weaker owing to the finite Rabi frequency of the microwave drive Ω.

Density-density correlations
We now consider the density-density correlations g (2) in the last resonator. The g (2)function is a quantity that describes the likelihood to measure two photons at the same place. The g (2) function of the last resonator is the normalized meanvalue of the second order moment of the field operators in the last resonator, Classical, thermal fields have g (2) -values larger or equal to unity with the coherent field exhibiting a g (2) -value of 1. A g (2) -value below 1, meaning that the photons are antibunched, is a sufficient condition to call the field quantum mechanical in the sense that there is no classical field showing the same results in measurements of the g (2) -function.
With recently developed refinements of microwave measurement techniques [31,32], measurements of g (2) -functions in circuit QED are now becoming feasible. In figure 6 we plotted the g (2) -function of the field in the last stripline resonator. To get a more detailed insight of the processes leading to a g (2) -value for specific parameters we plotted the g (2) -function along special values of the microwave drive frequency and the ratio E J /E C marked by white lines in figure 6. Figures 8 and 9 show the results for the different paths, denoted by a), b), and c) in the density plot of the g (2) -function in figure 6, for the g (2) -function as well as the corresponding values for the density of polaritons in the last cavity, c † +,2 c +,2 , and the second order moment, c † +,2 c † +,2 c +,2 c +,2 . For small values of E J /E C our system is basically linear because the nonlinearity U + /2 = (E C /2) sin 4 (θ) is negligible. A harmonic field mode driven by a coherent source is in a coherent state. Therefore the g (2) -function is equal to one for small values of E J /E C . As the nonlinearity grows for increasing values of E J /E C the g (2) -function plotted against the ratio of E J /E C and the frequency of the microwave drive ω µw gets more structured. In the density plot of the g (2) -function in figure 6 we can identify resonances where the frequency of the microwave drive matches the eigenenergies of the unperturbed system without microwave drive and dissipation. These resonances manifest themselves as separating lines between bunching regions (values of g (2) > 1) and anti-bunching regions (values of g (2) < 1).
To understand the origin of these separating lines, it is illustrating to analyze our system in terms of a symmetric mode, d + , and an antisymmetric mode, d − , where rather than the two localized modes c +,1 and c +,2 . In terms of d + and d − the Hamiltonian reads, The Hilbert space of the Hamiltonian H c+ can be described by two different bases, states that are labeled by the number of excitations in the collective modes, or states that are labeled by the number of excitations in the localized modes, The lines separating bunching and anti-bunching regions in figure 6 can now be identified with the energies of the 1 excitation states, Figure 7. Sketch of the energy spectrum of the Bose-Hubbard Hamiltonian H c+ for a two site model for vanishing nonlinearity U + compare a) and vanishing intersite coupling J + compare b). For vanishing nonlinearity a microwave drive can drive multiple transitions leading to a coherent state. Contrary to the linear case for strong nonlinearity one can only drive a transition between two distinct states as the energy differences between the eigenenergies aren't degenerate any more and the energy of a 2-excitation state, To understand the origin of the anti-bunching regions for a microwave drive that is blue detuned with respect to the energies of the states (13a-c) and the bunching regions for a red detuned microwave drive one has to consider the spectrum of the Hamiltonian H c + . For small nonlinearity, that is for values of E J /E C < 50, the Hamiltonian (12) reduces to a Hamiltonian for two uncoupled harmonic oscillators described by the modes d + and d − with energies ω + − J + and ω + + J + respectively. The eigenenergies in this situation are shown in figure 7 a). A microwave drive with frequency ω + −J + as depicted in figure 7 not only drives the transition from the groundstate to the first excited state of the symmetric collective mode |0 0 cm → |1 0 cm but also all other transitions to higher excited states |n 0 cm → |n + 1 0 cm . As a result the steady state in this situation is always the coherent state exhibiting a g (2) -value of 1. For slightly increased values of the nonlinearity that remain in the range U + < Γ c + , the system can still be described in terms of two weakly interacting collective modes. But the symmetric as well as the antisymmetric mode are subject to the nonlinearity and an intermode interaction, c.f. Hamiltonian (12). This can be seen in figure 8 a) where we plotted the g (2) -values that deviate from the value of a coherent field. The g (2) -function shows anti-bunching regions for blue detuned microwave drive with respect to the energies of the states (13a-c) and bunching regions for red detuned microwave drive. To gain insight into the underlying physical principles in this situation we calculated the density c † +,2 c +,2 and the second symmetric moment c † +,2 c † +,2 c +,2 c +,2 by an iterative meanfield approach to solve the master equation (11) with the Hamiltonian written as in (12). Operator mean values of a single driven dissipative mode with Kerr nonlinearity can be computed exactly [34] and we expand this model in a meanfield way to incorporate the denistydensity coupling. With this method we get good agreement with the numerical exact values for the density in the last cavity and are able to compute values for the polariton density close to the systems eigenenergies (13a-c) where our numerical approach fails to converge. For details about the method please see Appendix A. These results support our assertion that the system can be described by weakly interacting collective modes in the limit of small nonlinearities U + . In figure 8 a) numerically exact values are plotted in solid lines and values obtained by the above mentioned mean field method are plotted in dashed lines.
For strong nonlinearity U + and small intersite coupling J + , that is for values of E J /E C > 50, the c + -polaritons become transmon excitations and the Hamiltonian H c + splits into two parts describing the first and the second transmon qubit respectively. Here the collective modes d + and d − no longer decouple and the localized modes c 1 and c 2 become a more appropriate description of the system. The eigenenergy spectrum in this situation is shown in figure 7 b). The main difference to the spectrum without nonlinearity is that the microwave drive can't be adjusted to drive multiple transitions. In order to drive the transition to the state |02 s for example one has to adjust the microwave frequency to match half of the energy difference between the groundstate and the 2-excitation state |02 s because it is a two photon transition. Due to the anharmonicity of the eigenenergy spectrum no other transition can be driven. The difference of microwave frequencies needed to drive the transition from groundstate to |01 s respectively |02 s amounts to U + /2 which is bigger than the linewidth Γ c + . To get an estimate for the value of g (2) we simplify our model assuming that the frequency of the microwave drive is adjusted such that it resonantly drives a transition between the groundstate of our model |00 s and some excited state |0n s . Provided the Rabi frequency Ω and the loss rates κ and γ are all small compared to the frequency separation between different resonance lines, the system can then be modeled by a two level system consisting of the groundstate of our model |00 s and the excited state |0n s . In this situation the maximal occupation inversion one could get in the steady state is, ρ max = 1 2 (|00 00| + |0n 0n|) , and the g (2) -value for this density matrix would be g (2) = tr ρ max c † +,2 c † +,2 c +,2 c +,2 tr ρ max c † +,2 c +,2 2 = 2n(n − 1) n 2 which is below one for a 1-excitation state and above one for every state containing more than 2 excitations. Therefore bunching areas arise if states containing more than 2-exciations are excited and anti-bunching areas arise if only 1-exitation states can be excited and the photons pass the setup "one by one". In our Bose-Hubbard model the on-site nonlinearity is negative and hence all transition frequencies to states containing more than two excitations are red detuned with respect to transition frequencies to states containing only one excitation (13a-c). This is why bunching areas arise for red detuned microwave drive and anti bunching areas arise for blue detuned microwave drive. If we adjust the microwave drive frequency for every value of E J /E C to match the eigenfrequency of the antisymmetric 1-excitation state we get the transition from a perfectly uncorrelated field with g (2) = 1 to strongly correlated, anti-bunched field statistics with g (2) < 1 see figure 9. For a quantum phase transition of the ground state of the Bose-Hubbard Hamiltonian one would expect this transition as a consequence of the interplay of the intersite hopping J + and the on-site nonlinearity U + . For the driven dissipative system we observe that the interplay between the Rabi frequency of the microwave drive Ω cos(θ) and the on-site nonlinearity U + determines the particle statistics. This can be seen in figure 9 where we plotted the g (2) -function and the intersite coupling, on-site nonlinearity and the Rabi frequency of the microwave drive.

Summary
We have shown that a chain of capacitively coupled stripline resonators each coupled to a transmon qubit can be described by a Bose-Hubbard Hamiltonian for two species of polaritons. The validity of our approach has been checked for realistic parameters of the transmon qubits and stripline resonators and for low densities of polaritons. In a driven dissipative regime where a microwave source coherently drives the first cavity one can selectively excite only one species of polaritons and investigate the properties of a driven dissipative Bose-Hubbard model. We calculated the density and the g (2) -function of the polaritons in the last resonator of a two site setup and investigated their dependencies on the microwave drive, the intersite coupling and the on-site nonlinearity. For vanishing nonlinearity the g (2) -function is approximately equal to unity indicating a coherent field. With increasing nonlinearity bunching and anti-bunching areas arise depending on the frequency of the microwave drive. For a microwave drive that is in resonance with a transition to a state with one excitation, the polaritons are anti-bunched. If, on the other hand, the microwave drive can resonantly excite states containing more than 2 excitations in a multi-photon transition, the polaritons become bunched. If we adjust the microwave drive frequency to match one of the system's single excitation eigenenergies and compute the g (2) -function for different values of on-site nonlinearity,  Figure 8. Plots of the g (2) -function, the density of polaritons c † 2 c 2 and the second order moment c † 2 c † 2 c 2 c 2 in the second cavity for special values of E J /E C are shown. For all plots we have chosen the intersite coupling constant and the on-site nonlinearity to be J 0 /ω r = 0.04 and E C /ω r = 0.04 and the decay rates of transmon qubit and stripline resonator to be γ/ω r = 0.00008 and κ/ω r = 0.00004 respectively but we applied different Rabi frequencies of the microwave drive: for a) Ω/ω r = 0.004, and for b) Ω/ω r = 0.001. Plotted are results obtained by numerical calculation of the masterequation in solid lines and results obtained by a mean field approach with an exact single-site solution in dashed lines. Eigenenergies of the system without dissipation and driving are signalized by vertical dash-dotted lines. For E J /E C = 25 one can see clearly separated resonances for the symmetric and antisymmetric states d † ± |00 and a two photon resonance for the state d † + d † − |00 . The shape of the resonances at d † + |00 and d † − |00 are reproduced by the meanfield approximation and therefor a single mode feature. The two photon resonance for the sate d † + d † − |00 is not reproduced by the meanfield approach since it does not correctly incorporate the interactions between d † + and d † − modes. For E J /E C = 125 multiple resonances determined by the eigenenergies of the system without dissipation and driving arise .
intersite coupling and Rabi frequency of the microwave drive we see a transition from coherent to anti-bunched field statistics. That is, the polaritons are uncorrelated for small nonlinearity and exhibit a transition to anti-bunched behaviour as the on-site nonlinearity becomes larger than the Rabi frequency of the microwave drive. All our findings could be explored in experiments based on readily available technology.  Figure 9. A plot of the g (2) -function, the density of polaritons c † 2 c 2 and the second order moment c † 2 c † 2 c 2 c 2 in the second cavity is shown. The microwave drive frequency is chosen to match the transition from the ground state to the symmetric 1-excitation eigenstate. For all plots we have chosen the intersite coupling constant and the on-site nonlinearity to be J 0 /ω r = 0.04 and E C /ω r = 0.04, the decay rates for transmon qubit and stripline resonator to be γ/ω r = 0.00004 and κ/ω r = 0.00008 respectively and the Rabi frequency of the microwave drive to be Ω/ω r = Γ/ω r = 0.00004. Plotted are results obtained by numerical calculation of the masterequation in solid lines and results obtained by a mean field approach with an exact single-site solution in dashed lines. In resonance to the symmetric state g (2) shows a transition from uncorrelated coherent field particle statistics to anti-bunched correlated field particle statistics. The transition from coherent to anti-bunched is determined by the interplay of the Rabi frequency of the microwave drive and the on-site nonlinearity.