Photon condensation in circuit QED by engineered dissipation

We study photon condensation phenomena in a driven and dissipative array of superconducting microwave resonators. Specifically, we show that by using an appropriately designed coupling of microwave photons to superconducting qubits, an effective dissipative mechanism can be engineered, which scatters photons towards low-momentum states while conserving their number. This mimics a tunable coupling of bosons to a low temperature bath, and leads to the formation of a stationary photon condensate in the presence of losses and under continuous-driving conditions. Here we propose a realistic experimental setup to observe this effect in two or multiple coupled cavities, and study the characteristics of such an out-of-equilibrium condensate, which arise from the competition between pumping and dissipation processes.


Introduction
Massive particles and spin systems have been at the center of investigations in quantum many-body physics since the inception of the field. Photonic systems, however, have been the object of interest in this field only much more recently. This is, on the one hand, due to the fact that photons are intrinsically non-interacting, and on the other hand, due to the difficulty to confine and experiment with photons before they decay. One important exception in this context are ideas to produce a Bose-Einstein condensate (BEC) of weakly interacting polaritons [1,2,3,4,5], whose large photonic component is essential to achieve a low effective mass and therefore transition temperatures which are accessible in solid-state experiments. While current cold-atom experiments offer a well controlled setting to study bosonic quantum gases [6], the investigation of photonic condensates is still an active area of research, stimulated by new experimental directions in the field [7,8,9] and controversial discussions on the exact nature of such condensates [10,11,12].
More recently, the connection between many-body physics and photons has been addressed again from the perspective of photonic quantum simulators [13,14,15,16,17,18,19,20,21,22]. Here, the general goal is to implement strongly interacting many-body systems in a controllable way, to simulate the properties of non-trivial condensed-matter models. This can be achieved, for example, using ideas from cavity quantum electrodynamics (QED) [23], where effective photon-photon interactions can be obtained through the coupling to an intermediary system, such as atoms [24,25] or solidstate systems [26,27,28]. Based on this principle, various schemes for implementing bosonic Hubbard models for photons on a lattice [13,14], photonic quantum Hall systems [29,30,31,32], and strongly interacting photons in a 1D continuum [33,34], have been theoretically investigated. Although the experimental implementation of these ideas is still challenging, the development of scalable cavity-QED systems in onchip devices [35] is rapidly progressing, and analogous ('circuit QED') systems in the microwave regime [36,37,38,39,40,41] already approach the stage where photonic lattice models can be realized.
A central and still open question in the field of photonic quantum simulators is how to prepare and probe quantum many-body states of light. Because of unavoidable losses and the absence of a chemical potential (arising from an equilibrium particle reservoir), photonic many-body systems must necessarily be studied under driven and non-equilibrium conditions. Therefore, familiar concepts from condensed-matter physics, such as ground-state or equilibrium phase diagrams, are a priori not accessible in photonic many-body systems. In previous works it has been suggested to study, for example, the transient dynamics of an initially pumped system [15,17,42], or to use excitations spectroscopy and photon-correlation measurements of a weakly driven system to infer certain properties of the interacting photonic system [34,43,44,45,46]. In this work we consider a different scenario and study the dissipative dynamics of a strongly driven cavity array in the presence of engineered dissipation [47,48,49,50,51,52,53,54]. More precisely, we will show how in a circuit QED setting, the coupling of microwave photons to superconducting qubits can be used to scatter photons from high to low momentum states, eventually accumulating photons in a zero-momentum condensate. Previously, a similar scheme has been explored as a dissipative way to prepare a Bose Einstein condensate (BEC) of atoms [48]. In the case of photonics systems, it is essential that this mechanism implements a numberconserving coupling of photons to an effective low-temperature bath, and therefore provides a new approach for the preparation of quasi-equilibrium many-body states in open and driven photonic systems. Here we propose and analyze a proof-of-principle experiment to study condensation of microwave photons in a dissipative cavity array, and describe the properties of such a non-equilibrium BEC, which arise from the interplay between driving, decay, and thermalization processes. Such experiment would realize a controlled setting for the simulation of non-equilibrium condensation phenomena, and more generally, would offer a new route towards preparing stationary states of strongly interacting photons. Moreover, our work shows that the flexibility provided by circuit QED, which so far has mainly been employed to engineer strong coherent interactions, can equally well be used for the design of various non-trivial dissipative couplings, and be applied to simulate driven quantum systems in unconventional environments. Fig. 1 illustrates the basic setup for an array of L coupled cavities, where each cavity is represented by a single photonic mode of frequency ω c and bosonic annihilation operator c ℓ . Photons can tunnel between neighboring cavities with hopping amplitude J and, in addition, local interactions with two-or many-level systems can be used to induce Kerrtype nonlinearities with an effective photon-photon interaction strength U. A generic model for this system is then given by the Bose-Hubbard Hamiltonian (see e.g. Ref. [13])

Coupled cavity arrays
Various generalizations of this model have been discussed in the literature and several implementations have been proposed using optical or microwave cavities. The photonic Hubbard model given in Eq. (1) describes a gas of interacting bosons on a lattice, where the hopping J competes with on-site interaction U. However, in contrast to the coldatom physics scenario, the bosons here can decay, and under experimentally-relevant conditions ω c ≫ J, U, k B T (where T is the temperature and k B the Boltzmann constant), the equilibrium state of this model is simply the vacuum state. Therefore, in photonic many-body systems we are mainly interested in the out-of-equilibrium dynamics ofĤ c in the presence of losses and external driving fields. In particular, in this work we model the resulting dissipative dynamics for the system density operator ρ by a master equation (ME) of the form, where ℓ ) describes an external driving field of frequency ω d which is used to excite the system, and the second term accounts for photon losses in each cavity with a field decay rate Γ. While a finite driving field is required to counteract the losses, it will in general also compete withĤ c and, for strong driving fields, even dominate the system dynamics. Therefore, in previous works it has been suggested to study either the transient dynamics of an initially prepared photonic state [15,17,42] (where Ω ℓ = 0 for times t > 0) [42] or to use weak excitation spectroscopy [34,43,44,45,46] (Ω ℓ → 0) to probe the many-body spectrum of HamiltonianĤ c .
In this work we are interested in the opposite regime of a strongly and continuously driven cavity array. We study the dynamics of this system in the presence of an additional artificial thermalization mechanism, denoted by L κ in Eq. (2). More precisely, we will show below how a non-local coupling of photons to superconducting qubits can be engineered in an array of microwave cavities to implement a dissipative photon scattering process of the form The interpretation of this term can be seen best in the case of just two cavities. Then, the first term in Eq. (3) describes the scattering of photons from the asymmetric (energetically higher) modeĉ a ≡ (ĉ 1 −ĉ 2 )/ √ 2 into the symmetric (energetically lower) modeĉ s ≡ (ĉ 1 +ĉ 2 )/ √ 2, while preserving the total photon number. The second term in Eq. (3) describes the reverse process. In the absence of losses, the action of L κ would, for two cavities, drive arbitrary initial photon states into a thermal equilibrium distribution, characterized by a detailed balance between symmetric and antisymmetric modes, with the ratio κ ′ /κ = exp(−2J/k B T eff ) defining an effective temperature T eff . Similarly, for a whole cavity array, the processes ∼ κ scatter photons to lower momenta and -as previously shown for an analogous system of cold atoms [48] -drive the system towards a condensate of photons in the zero momentum mode, |BEC = (ĉ † q=0 ) N |0 / √ N!, where q is the quasimomentum and N the total number of particles. Also in this case the opposite processes ∼ κ ′ can be roughly interpreted as a finite temperature effect, although, as it will be shown below, for L > 2 and κ ′ = 0 the stationary state of L κ does no longer represent a true thermal equilibrium state.
The competition of L κ with a Kerr-type nonlinearity in the Hamiltonian has already been discussed in the atomic, number-conserving setup [55,56]. However, in the photonic case, under realistic conditions, L κ competes with external driving fields and intrinsic photon losses as described by the full ME (2). To study this competition more explicitly, we focus in the following on the limit U → 0, and consider the experimental scenario where neighboring cavities are driven by coherent fields of alternating amplitude Ω ℓ = (−1) ℓ Ω (see Fig. 1). In this case, the accumulation of photons in the q = 0 mode can be interpreted as a clear signature for photon condensation induced by the dissipative photon-photon interactions in L κ .

Physical implementation
In this section we show how the photon scattering mechanism described by Eq. (3) can be implemented in an array of superconducting microwave resonators. The basic setting is illustrated in Fig. 2a for two cavities, each coupled to a nonlinear element, for example a superconducting Cooper pair box ('charge qubit'). The qubits are placed next to each other to obtain strong electrostatic or magnetic interactions. As we discuss now, this configuration allows us to engineer both nonlinear as well as non-local dissipation processes for photons.
For simplicity we restrict the following analysis to a single block of only two cavities, but the generalization to a whole array of linked cavities is straightforward. The Hamiltonian for this system can be derived from the corresponding equivalent circuit model schematically shown in Fig. 2a. By restricting each resonator to a single modê c ℓ , ℓ = 1, 2, and by approximating each Cooper pair box by a two-level system with ground state |g and excited state |e , we obtain whereσ (ℓ) − ≡ |g ℓ e ℓ |. The first part is the free cavity Hamiltonian, while the second term describes the dipole interaction between photons and qubits with strength g. Finally, . Two stripline resonators are coupled through a system of two mutually interacting Cooper pair boxes. b) Four-level diagram corresponding to the two coupled qubits, whose dynamics is described by Eq. (6).
with rate 2Γ and qubit decay with rate γ. The full system dynamics is then described by the MEρ and a summary of the energy levels and the most relevant transitions is shown in Fig. 2b. Our goal in the following is to eliminate the qubit dynamics and derive an effective ME for the cavity modes. Specifically, we are interested in the dissipative two-photon process, where the qubits change from |A to |S by absorbing a photon from the antisymmetric cavity mode (ĉ a ) and emitting a photon into the symmetric cavity mode (ĉ † s ) (see Fig. 2b). Since the overall process is O(g 4 ) in the photon-qubit coupling strength, a general derivation is quite involved, and to make the following discussion more transparent we will now focus explicitly on the hierarchy of energy scales drawn in Fig. 2b. In particular, we assume that δ, ∆ e , and ∆ ≡ δ − ∆ e are much larger than all the other frequency scales g, Ω a , ∆ s,a , Γ, γ. In this limit, none of the single-photon processes is resonant and we can use a Schrieffer-Wolff transformation to derive an effective two-photon Hamiltonian. We perform the unitary transformationH ≡VĤV † withV = eŜ, and make the ansatẑ We defineĤ 0 andĤ g as the first and last lines of the Hamiltonian (6), respectively, and choose the coefficients α i such that [Ŝ,Ĥ 0 ] = −Ĥ g . This can be achieved by setting In view of the assumptions discussed above, we will use the approximate results α 1,η ≈ g/δ, α 2,η ≈ g/∆ and define the parameter ǫ = ∆/δ < 1. After this transformation the Hamiltonian (6) reads HereH q is the qubit HamiltonianĤ q written above, with small frequency shifts absorbed into a redefinition of the detunings ∆ e,s,a . The modified cavity Hamiltonian is and the average · is taken with respect to the stationary qubit state in the limit g → 0. For the parameter regime of interest, P sa ≈ 0, andH c corresponds to the free cavity Hamiltonian with a qubit-mediated tunneling amplitude J. Finally, the effective coupling between photons and qubits can be written as Note that inH int we have only kept resonant two-photon processes and already omitted terms like ∼ gΩ a /∆ ×ĉ s,a and g 2 /∆ ×ĉ 2 s,a . Both of these terms oscillate with the detuning δ and will only give small corrections to the results presented below.
The transformationV induces a weak mixing between photon and qubit states, which apart from generating new effective interactions also modifies the dissipative couplings. Since we consider γ ≫ Γ, we find that the main result of this mixing is an enhancement of the cavity decay rates Γ → Γ a,s ≡ Γ +Γ s,a . After a transformation of the jump operators in the ME (7), and averaging over the qubit states, we find that for the parameter regime of interest, these rates arẽ In summary, we find that the system dynamics in the new dressed state representation is well described by the MĖ where In the limit g eff → 0, qubits and cavities are decoupled, and the system relaxes on a timescale γ −1 into the state ρ(t) ≃ ρ c (t) ⊗ ρ 0 q , where L q ρ 0 q = 0. For finite g eff < γ, we can use a perturbative projection operator technique to eliminate the qubit degrees of freedom and derive an effective ME for the reduced cavity density operator [59], Evaluating this expression we obtaiṅ wheren η =ĉ † ηĉ η andĤ eff is an effective photon-photon interaction In Eq. (22) we have defined κ ≡ g 2 eff Re{S as,sa (0)}, Γ Φ ηη ′ ≡ g 2 eff Re{S ηη,η ′ η ′ (0)}, and the interactions U eff ≡ g 2 eff Im{S as,sa (0)} and U Φ ηη ′ ≡ g 2 eff Im{S ηη,η ′ η ′ (0)}, given in terms of the qubit correlation spectra which can be calculated using the quantum regression theorem [60].

Discussion
By ignoring coherent photon-photon interactions for the moment, we see that the first line of the effective cavity dynamics given in Eq. (22) represents the dissipatively coupled cavities as introduced in Eq. (2). In particular, we can write By comparison of Eqs. (22) and (3) we identify κ ′ = ǫ 2 κ, with the parameter ǫ 2 ≡ (∆/δ) 2 determining the minimal effective temperature of the engineered two-photon process. For typical parameters, ω q = 10 GHz, δ = 1 GHz and ∆ e = 0.8 GHz, we obtain ǫ = 0.2, and hence the limit κ ′ → 0 considered in most parts of the paper can be implemented to a very good approximation. However, we point out that the ratio κ ′ /κ can always be increased artificially, for example, by adding an additional coherent or incoherent driving field to populate the state |S . Our derivation also shows the existence of coherences betweenĉ † sĉ a andĉ † aĉ s , which lead to additional squeezing terms of order O(κǫ), and that in Eq. (25) are summarized by L ǫ . For ǫ ≪ 1, and in the presence of cavity losses, we expect the influence of these coherences to be negligible, while the relevant effect on the populations is already captured by a finite κ ′ . Indeed, numerical simulations similar to those discussed in Sec. 4 show no significant differences between the exact ME (22) and our model Liouvillian (3), even up to values of ǫ ≈ 0.5. However, by choosing ǫ ≈ 1 and Γ ≪ κ, it would be possible to extend our model and study physical effects directly linked to the parameter ǫ, such as number-conserving squeezing as previously analyzed in [61].
Apart from the two-photon scattering processes of interest, the coupling to the qubits introduces various imperfections. On the one hand, these are the enhanced cavity decay rates Γ s,a introduced above, and on the other hand, we obtain additional dephasing terms Γ Φ ij due to fluctuating Stark shifts when the qubit jumps incoherently between states |G and |A . Since the ratio κ/Γ s,a scales as ∼ g 2 /γ 2 , single photon losses can be suppressed in the strong coupling regime, as long as κ also exceeds the bare decay rate Γ. The ratio κ/Γ Φ ij is independent of g, and depends on the details , as a function of qubit driving Ω a and detuning ∆ a . Rest of parameters as in (a). c) The same plot as in (b) but for low-momentum modes, of the correlation functions (24) and therefore on the parameters Ω a , ∆ s,a . In Fig. 3a, κ is compared with the other decoherence rates, and we see that an experimentally feasible parameter range exists, where κ is the dominant rate. The main corrections arise from the enhanced decay Γ a and the dephasing Γ Φ aa of the antisymmetric mode, while the corresponding rates for the symmetric mode are much smaller. In Fig. 3b and 3c we plot the two quantities F η ≡ κ/(Γ η + Γ Φ ηη ), η = s, a, as a function of the control parameters Ω a , ∆ a . These quantities represent a simple figure of merit for characterizing the validity of our model for low (F s ) and high (F a ) momentum modes. We see that within a large parameter regime, the number conserving photon scattering rate κ can exceed all decoherence processes discussed here, and therefore the proposed model in Eq. (2) provides a reliable description of the system dynamics. In particular, this is true for the low momentum (symmetric) modes, where condensation occurs.
Finally, let us briefly comment on the effective photon-photon interactions described byĤ eff . For a resonant two-photon process, ∆ s = 0, we find U eff ≈ 0. Moreover, U Φ ηη ′ ≈ Γ Φ ηη ′ < κ for the typical parameters considered above, and in particular these interactions vanish for ∆ a = 0. Therefore, in this work we focus exclusively on the purely dissipative cavity dynamics. However, by setting |∆ s | > γ instead, the regime of strong coherent photon-photon coupling U eff > κ > Γ Φ ηη ′ , Γ s,a can be engineered as well, havingĤ eff ≃ U effĉ † sĉ sĉ † aĉ a .
In summary, we find that the ME (2) can be implemented with superconducting microwave cavities under realistic experimental conditions. Our analysis, which has been detailed here for generic coupled two level systems, can be thus adapted to various nonlinear Josephson devices, where the parameter λ x,z plays the role of the nonlinearity.

Semiclassical approximation
In the previous section we have shown how to implement the effective ME (2). Now, we discuss how the scattering between photons induced by the Liouvillian L κ leads to condensation of photons in an open and driven cavity array. In this section we first illustrate the basic process for the minimal case with L = 2 cavities only. As already mentioned above, we assume U = 0, and consider the case where only the antisymmetric mode is excited, i.e. Ω 1 = −Ω 2 = Ω/ √ 2. In the rotating frame with respect to the drive frequency ω d the equations of motion (EOM) for the field expectation values ĉ s and ĉ a of the symmetric and antisymmetric cavity modes are where δ d = ω c + J − ω d is the detuning of the antisymmetric mode with respect to the driving field. For the moment, we assume κ ′ = 0 in Eq. (3) and postpone the discussion of this term to the end of the present section. The EOM for the field expectation values couple to third-order correlation functions, starting a hierarchy of equations that cannot be solved analytically. To truncate the hierarchy, we resort here to a semiclassical approximation in which the state of the system is assumed to be coherent. With this assumption, each normal-ordered correlation function can be readily substituted by the corresponding product of field expectation values. The EOM read where s ≡ ĉ s and a ≡ ĉ a for brevity. The semiclassical approximation, in general, is valid for large photon numbers but it offers a satisfactory explanation of the dynamics of the system for smaller occupation as well, as we verify below with the numerical integration of the EOM (26). The exact equations (26) become nonlinear c-number equations for the expectation values s and a. Nonlinear equations may have multiple, stable and unstable steady-state solutions. A paradigmatic example of this behavior are the classical equations of the laser [62], where a solution which lacks coherent emission is possible at any pumping strength, but is unstable above a certain threshold. Here we find a comparable result, with a threshold value for the pumping strength Ω. Below the threshold, the only stable steady-state solution to Eq. (27) is (c.f. also Appendix A) and only the pumped antisymmetric mode is populated. For driving strengths above the threshold we obtain the stable solution where the symmetric mode features nonvanishing occupation. Note that these equations depend on J only via the detuning δ d , which has the only effect to reduce the driving strength. Therefore, from now on we will restrict ourselves to δ d = 0 and J → 0. The stationary solution to the semiclassical EOM then simplifies to The behavior of the populations n s ≡ |s| 2 and n a ≡ |a| 2 is plotted in Fig. 4 as a function of the pumping intensity and the effective decay κ. Below threshold, the symmetric mode is empty and the antisymmetric mode occupation grows quadratically with the pumping intensity. For Ω > Ω crit , the antisymmetric mode population saturates, and the symmetric mode occupation grows linearly. Thus, if Ω is sufficiently large, most of the photon population is transferred to the symmetric mode. This transition to a state in which photons 'macroscopically' occupy the symmetric superposition of two cavity modes is a direct consequence of the engineered nonlinear Liouvillian (3). The discussion of Eqs. (31) illustrates both, similarities and differences between the present setup and a laser [63,64,65,66]. In a single-or multi-mode laser, the light field is coupled to a pumped reservoir and the lasing threshold corresponds to the point where the total light intensity starts to grow. In our case, since the antisymmetric mode is driven directly by a classical field, the total photon population is always finite and grows monotonically. The threshold then marks the point where the many-body scattering mechanism contained in L κ dominates, and a macroscopic transfer of photons from the antisymmetric to the symmetric mode occurs. In this sense, one can rather think of the antisymmetric mode as a photon reservoir from which a condensation of photons ('lasing') into the symmetric mode occurs. However, as we will see below, this analogy between lasing and condensation [10,11,12,67,68] holds in our system only for certain parameter regimes, and in particular striking differences can be observed in the regime of low photon numbers.
Finally, note that both for κ → 0 as well as for κ → ∞ the condensation is suppressed, and that the phase diagram in Fig. 4 is symmetric with respect to the ratio κ/Γ. For small κ, not enough photons are scattered into the symmetric mode, while for large κ the antisymmetric mode is highly overdamped and the population of the whole system is small as in the first place.

Exact diagonalization
To support the results of the semiclassical approximation discussed above, we resort to the exact diagonalization of the system composed by two cavities. The time-evolution of the system for typical parameters features an initial transient in which the photonic population accumulates in the antisymmetric mode, which is directly pumped, and a later stage in which photons scatter into the symmetric mode. Finally, a steady state is reached in which a dynamical equilibrium takes place between the populations in the two modes. We remark that the population in the modes is continuously subjected to losses.
The steady-state values of the populations n s = ĉ † sĉ s and n a = ĉ † aĉ a in the symmetric and antisymmetric mode, respectively, are shown in Fig. 5 as a function of the pumping strength Ω and for two different values of κ/Γ. In the limit κ ≪ Γ, where we expect the transition to occur at large mode occupation numbers n a , we find that the exact solution matches qualitatively and quantitatively very well the semiclassical predictions. In particular, we see that above Ω crit the population n s grows linearly with the driving strength, while the antisymmetric mode saturates at a value n a = Γ/κ. Differently from the prediction of the semiclassical analysis, however, n s does not vanish exactly below the threshold but it is smaller than n a . A study of the equal-time twophoton correlation g (2) ≡ n s (n s −1) / n s 2 reveals that, when crossing Ω crit , the photon statistics of the symmetric mode changes from a thermal (g (2) ≈ 2) to a Poissonian (g (2) ≈ 1) distribution. This aspect is in agreement with the standard lasing transition.
In the opposite regime, κ ≫ Γ, Eq. (31) predicts that the population n a of the driven mode is always less than one, and we do not expect the semiclassical analysis to give accurate results. Indeed, Fig. 5b shows that in this limit the transition is completely washed out and n s exceeds n a for all pumping strengths. Nevertheless, above Ω crit we still observe the characteristic saturation of n a and -apart from an additional offsetthe correct linear scaling of the population in the symmetric mode. However, in strong contrast to the semiclassical regime, we find that the photon statistics of the symmetric mode is sub-Poissonian (g (2) < 1) at low driving strengths and approaches the classical limit g (2) = 1 from below. This anti-bunching effect can be understood from the fact that the effective damping of theĉ a mode is given by Γ + κ( n s + 1). Therefore, the scattering of the first photon into the symmetric mode changes the damping significantly and suppresses the following repopulation of the 'reservoir mode'ĉ a . The limit κ ≫ Γ and Ω ≈ Ω crit would then result in a 'single photon condensate' with n a ≪ n s ∼ 1 and a non-classical statistics with g (2) → 0.5.

Effective temperature
Finally, let us briefly comment on the effect of the reverse photon scattering process ∼ κ ′ in Eq. (3). As we have discussed in Sec. 2, in the case of two isolated cavities the presence of the reverse photon scattering term can be interpreted as a finite temperature effect, and we can set κ ′ = κe −2J/k B T eff . By including this term, the EOM are These equations are of the same form of (26), with the replacements Γ → Γ + κ ′ and κ → κ − κ ′ . Therefore, in the semiclassical regime we obtain the same physics, but with a renormalized pumping threshold In particular, for any fixed driving strength we can use Eq. (33) to define a critical effective temperature above which photon condensation does not occur.

Photon condensation in a cavity array
We turn now to the solution of Eq. (2) in a lattice with L ≫ 1 cavities. The central goal of this section is to show that the many-body Liouvillian (3) induces condensation of photons into the homogeneous lattice mode with zero momentum q. Following the scheme introduced in the previous section in the case of L = 2 cavities only, we consider the scenario in which each cavity modeĉ ℓ in the array, with ℓ = 1, . . . L, is coherently driven by a classical field with staggered amplitude Ω ℓ = (−1) ℓ Ω. In this way, the coherent driving acts on the edge of the Brillouin zone, while dissipation-induced condensation is expected to take place at the center. For simplicity, we assume periodic boundary conditions for the lattice. We rewrite Eq. (2) in terms of the annihilation operatorsĉ q = 1 √ L ℓ e iqℓĉ ℓ of the photonic modes in momentum space, where q takes values −π(L − 2)/L, . . . , π(L − 2)/L, π (given here in units of the inverse lattice spacing) in the discretized Brillouin zone. In the rotating frame with respect to the coherent driving, the Hamiltonian is then given bŷ where Ω q = 1 √ L ℓ e iqℓ Ω ℓ and δ d = ω c + 2J − ω d is the detuning of the q = π mode from the driving frequency. The Liouville operators are and In the latter equation, the function K plays the role of a scattering kernel, which has the general form In the following analysis, we will focus on the case κ ′ = 0. Furthermore, as discussed above, the photon condensation effect does not depend on the tunneling amplitude and we can assume J → 0 in the Hamiltonian. Solving the ME (2) in an extended array is in general a very demanding task, and exact diagonalization can only be provided for very limited system size (see our results in Fig. 6b for L = 4 cavities). For this reason, we resort in the following to two complementary approximations, which allow us to treat the full dynamics of the photonic population in the whole Brillouin zone. First we consider the equation of motion for the field expectation values within the semiclassical approximation, which was discussed in Sec. 4 for the case L = 2, and showed reliable results for large photon numbers. Then we treat the photonic population directly and derive a Boltzmann-like equation starting from the exact EOM for the two-point correlation functions. Within the latter approximation we are able to estimate the finite correlation length of the photonic field in the array. In this respect, this technique fills the gap between the semiclassical analysis of an infinitely extended array and the numerical solution to a small-size system. We emphasize that all techniques agree on the central results that we present here, i.e. the existence of a macroscopic photonic population at the center of the Brillouin zone.

Semiclassical approximation
For many cavities, the EOM for the cavity fields ψ q = ĉ q (analogous to Eq. (27)) read where Ω q = δ q,π Ω π . The numerical solution of this system of nonlinear coupled equations is straightforward. Fig. 6a shows the population n q = |ψ q | 2 on the discretized Brillouin zone for L = 22 cavities. We see that, in the initial stage of the dynamics, the coherently pumped mode at q = π is populated on the time scale of Γ −1 . Later, photonic population is transferred to the mode q = 0 and the steady state is reached. A surprising feature is that the population of the intermediate modes in the Brillouin zone is almost vanishing (and indeed not visible in Fig. 6a) throughout the time-evolution. This contrasts sharply with typical relaxation phenomena in the Brillouin zone (see e.g. Ref. [69]) where the excited population is expected to drift to lower energy states by interaction-induced relaxation. The reason for this behavior is the specific form of the engineered systembath coupling L κ . We can take advantage of the negligible population in the intermediate modes of the Brillouin zone and introduce a two-mode approximation, in which only the averages ψ 0 and ψ π of the modes q = 0 and q = π, respectively, are assumed finite. The EOM (38) then reduce to The structure of these equations is exactly the same as Eq. (27), with the |a| 2 /2 and |s| 2 /2 replaced by the photon densities |ψ π | 2 /L and |ψ 0 | 2 /L, respectively. The analytical discussion of these coupled equations follows the same path as that from Eq. (27) to Eq. (31). The threshold, above which photons start to condense into the q = 0 mode, is Ω π,crit = L/2 × Γ(Γ + 2κ) 2 /2κ. Note that compared to the results presented in the previous section, a factor 2κ, instead of κ, appears. This is simply due to the fact that in a 1D array (with periodic boundary conditions) each cavity is dissipatively coupled to two neighboring sites.
In Fig. 6b we also present the results of the exact diagonalization of a small cluster of L = 4 cavities, for typical parameters comparable to the semiclassical results. We show the absolute values of the entries of the one-particle density matrix in momentum space, ĉ † qĉ q ′ . The values along the diagonal are the populations n q , and the peak at q = 0 clearly demonstrates substantial occupation of that mode. The non-diagonal entries, which represent the correlation between modes with different momentum, are small. These values cannot be obtained correctly from the semiclassical approximation, where correlations factorize and one would have | ĉ † qĉ q ′ | ≃ √ n q √ n q ′ . The absence of correlation between modes with different momentum can be understood, since the photonic scattering term originates from the non-unitary contribution to the ME (2). This fact suggests that an approximation scheme should hold, where the correlations between momentum modes are neglected from the beginning, and only the populations n q are the dynamical variables. The Boltzmann-like EOM that implements this scheme is presented below.

Steady-state photon distribution
The exact EOM for the populations n q = ĉ † qĉ q are where the term S(q) contains fourth-order correlation functions, which are not diagonal in momentum space, and reads One route to reduce the EOM to a manageable form is to implement a truncation in the infinite hierarchy of coupled correlation functions, expressing higher-order correlations in terms of products of lower-order ones. Here, we reduce the fourth-order correlation function to a product of second-order diagonal correlation functions according to the prescription This prescription is in fact a Hartree-Fock approximation [70], because we assume that the dominant contribution to the correlations arises in the density channel of the momentum modes. The motivations for this choice have been discussed above. Using this factorization we obtain the following simplified form for the scattering term A self-consistent solution is to assume that ĉ q = 0 for q = π. The EOM (40) become then a set of coupled, nonlinear, Boltzmann-like equations [69] for the populations only. However, on the edge of the Brillouin zone, the coherent pumping Ω π couples to the average field, which does not vanish. The numerical integration of the resulting equations is straightforward and produces the distribution n q in the whole Brillouin zone. A finer discretization in the Brillouin zone is possible, with respect to the semiclassical equations (38), with a comparable numerical effort. Typical profiles of the steady-state distributions are shown in Fig. 7. Due to the increased resolution in momentum space, the broadening of the condensate peak at q = 0 can be appreciated, and its dependence on the parameters is discussed below.
To discuss the width of the central peak, which is related to the correlation length ξ of the condensate (which we give here in units of the lattice spacing), we can restrict our analysis to the dynamics of the π mode and to the populations of modes with q ≈ 0. For the coherence and population of the q = π mode we obtain and Here we have defined with ′ q excluding q = π, and used the Hartree-Fock decoupling approximation, which for three operators reads ĉ † q 1 +q 2 −qĉq1ĉq2 ≃ n q 1 δ(q 2 −q) ĉ q +n q 2 δ(q 1 −q) ĉ q . For δ d = 0, the normalized steady-state population N π = ĉ † πĉ π /L of the driven mode is then given by For the low-momentum modes we find and the steady-state solution can be approximately written in the form where Equation (49) is a Lorentzian profile with broadening 1/ξ. This scale plays the role of a correlation length since, in the translationally-invariant case, Here we have used that, in momentum space, the single-particle density matrix is mainly occupied on the diagonal. This approximation is in line with the Hartree-Fock approximation used above and can be physically justified by the fact that lowmomentum modes are populated by photons which are scattered incoherently from the q = π mode.

Discussion
For the Lorentzian low-momentum distribution (49) we obtain N 0 ≃ n 0 /ξ, and we can interpret the quantity L × N 0 as the total number of condensed photons, while the population L × N π serves as finite photon reservoir. In steady state, the values of N π , n 0 and ξ can be determined from the nonlinear Eqs. (47) and (50). The stability of the low-momentum modes requires that N π < Γ/4κ, but for a driving strength above Ω π,crit the value of N π will approach this limit. Therefore, in this regime whereΩ π = |Ω π |/Ω π,crit > 1. We see that above the critical driving strength the number of photons in the low-momentum modes increases linearly withΩ π . At the same time the correlation length increases ∼Ω 2 π , indicating an even stronger tendency for photons to occupy the q = 0 mode. Note that while in the classical high-photon limit, Γ ≫ κ, we obtain the ratio N 0 /N π ≃ (Ω π − 1), in the opposite case, κ ≫ Γ, we get an additional enhancement factor, N 0 /N π ≃ (2κ/Γ)(Ω π − 1). Therefore, in this regime, very pure, dilute photon condensates with N 0 ≈ 1 and ξ ≃ (κ/4Γ)(Ω π −1)Ω π ≫ 1 can be prepared. For typical parameters, the correlation length ranges from a few to a few tens of cavities, and can easily be tuned by adjusting the driving strength. This means that with realistic settings of around ten cavities, the variation of the correlation length could be studied in experiments.
Finally, we find it worthwhile to compare the steady-state distribution of the nonequilibrium open system studied here to a thermal equilibrium gas of massive bosons with finite chemical potential. To do so, we compare the distribution (49) with a Bose-Einstein distribution n BE q = 1/(exp{(ε q − µ)/k B T } − 1) (c.f. Fig. 7b), where ε q = −2J cos(q) in our case. Considering small momenta, a low q expansion of both distributions gives the following effective temperature and chemical potential for our photon condensate Note that the temperature T defined in this way is different from T eff associated to a finite κ ′ , and arises purely from the competition between photon losses and equilibration processes in a driven system. Above threshold, where n 0 , ξ ≫ 1, we obtain the scaling k B T /J ≃ 4/Ω and k B T /J ≃ 2Γ/(κΩ) for the limits Γ ≫ κ and Γ ≪ κ, respectively. This confirms the conclusion from above that, for κ/Γ → ∞, arbitrarily pure photonic condensates can be prepared with our scheme.
As a final remark, we point out that the detuning δ d of the driving field, which enters as an energy offset in the cavity Hamiltonian (34), affects the effective chemical potential µ only indirectly via a modification of Ω π,crit . In particular, and in strong contrast to equilibrium systems, the sign of δ d does not play a role in our non-equilibrium scenario and the chemical potential is related to the strength, rather than to the detuning of the driving field, as seen from the combined Eqs. (53) and (52). While for the present system this relation can be derived explicitly, it also suggests a general way to think about the effect of driving fields on the stationary states of open photonic quantum many-body systems.

Conclusions and outlook
In conclusion, we have analyzed a scheme to achieve condensation of microwave photons in an open and driven array of superconducting cavities. In particular, we have shown how the coupling to superconducting qubits can be used to engineer non-local and number-conserving dissipation processes for microwave photons, which mimic the coupling to an effective low-temperature bath. Under state-of-the-art experimental conditions, these processes can exceed the intrinsic losses in the system and produce stationary many-body states of photons, which are independent of the details of the external driving fields. We have proposed a basic experimental setup for demonstrating a stationary condensate of photons in two or multiple coupled cavities and compared the properties of such an out-of-equilibrium photon condensate to an equilibrium scenario for massive bosons.
More generally, our work illustrates the potential of superconducting cavity arrays for simulating various dissipative and out-of-equilibrium quantum many-body problems. The present example of photon condensation already shows how in engineered and fully controlled systems, an appropriately designed coupling to a bath can lead to nontrivial stationary states of an otherwise non-interacting photonic system, and could provide new insights for related experiments in condensed-matter systems. Beyond this basic example, the proposed scheme for engineering dissipation can be easily extended and combined with strong coherent photon-photon interactions. Compared to original ideas described in the context of cold atoms, superconducting circuits represent a complementary platform with, in many respects, additional flexibility to design nonstandard dissipation processes, as well as to study open, non-equilibrium quantum many-body systems. Figure A1. a) Fraction of photons in the symmetric mode, n s /n ≡ n s /n as a function of the dissipative rate κ, for incoherent pumping p/Γ = 0.5. b) Threshold given by Eq. (28), below which (29) is the only stable solution and above which (30) is the only stable solution.
Solving for the steady-state solution of these equations within the semiclassical approximation n sna ≈ n s n a , we find the condition p Γ = n n+2 , being n ≡ n s + n a the total number of particles in the steady state. Substituting then n a = n − n s = 2p Γ−p − n s in the first steady-state equation, for p = Γ, we find a smooth crossover between an equal population of symmetric and antisymmetric modes for κ ≪ p, Γ, and condensation to the symmetric mode for κ ≫ p, Γ (see Fig. A1a). Interestingly, for p → Γ, however, all photons occupy the symmetric mode for all values of κ = 0. In the latter case, heating (incoherent pumping) compensates cooling, and only the engineered dissipative mechanism survives.
Stability of the steady-state solutions. The stability of the different solutions can be checked by linearization of the EOM around the steady-state values (29) and (30). To this end we linearize the EOM (27) as s → s 0 + δs, a → a 0 + δa (similarly for the conjugate fields), being s 0 , a 0 the steady-state solutions and δs, δa representing fluctuations. This gives d dt (δs, δa, δs * δa * ) T = J(δs, δa, δs * δa * ) T with the 4 × 4 Jacobian matrix J: , whose eigenvalues will have negative real part in the region of stability. Evaluating the eigenvalues of this matrix at the steady state solutions, we find that the steadystate solutions presented in Section 4 are stable in the whole range of parameters. In Fig. A1b we show the threshold Ω crit given by Eq. (28), below which (29) is the only stable solution and above which (30) is the only stable solution.