You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: '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.'
address:
- '$^1$Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, A-6020 Innsbruck, Austria'
- '$^2$Institute for Theoretical Physics, University of Innsbruck, A-6020 Innsbruck, Austria'
author:
- 'D Marcos$^1$, A Tomadin$^1$, S Diehl$^{1,2}$ and P Rabl$^1$'
title: Photon condensation in circuit QED by engineered dissipation
---

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 [@DHY2010; @DWSBY2002; @KRKBJKMSASSLDD2006; @BHSPW2007; @ASLBVMLBKSTV2009], 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 [@BDZ2008], the investigation of photonic condensates is still an active area of research, stimulated by new experimental directions in the field [@KVM2010; @KSVM2010; @KSDVW2011] and controversial discussions on the exact nature of such condensates [@Snoke2003; @SL2002; @SLS2003].

More recently, the connection between many-body physics and photons has been addressed again from the perspective of photonic quantum simulators [@HBP2008; @TF2010; @HBP2006; @GTCH2006; @ASB2007; @RF2007; @NUTY2008; @FOU2008; @KH2009; @SGHBT2010]. 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) [@RBH2001], where effective photon-photon interactions can be obtained through the coupling to an intermediary system, such as atoms [@I1997; @BBMBNK2005] or solid-state systems [@HBWGAGFHI2007; @FFESPV2008; @RTM2009]. Based on this principle, various schemes for implementing bosonic Hubbard models for photons on a lattice [@HBP2008; @TF2010], photonic quantum Hall systems [@CAB2008; @Koch2010; @Hafezi2011; @UC2011], and strongly interacting photons in a 1D continuum [@CGMVLD2008; @CGTDCI2009], have been theoretically investigated. Although the experimental implementation of these ideas is still challenging, the development of scalable cavity-QED systems in on-chip devices [@LTHPH2011] is rapidly progressing, and analogous (‘circuit QED’) systems in the microwave regime [@WSBFHMKGS2004; @CBSNHM2004; @BHWGS2004; @SG2008; @CW2008; @MWBLLNCSWWYYZMC2011] 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 [@HBP2006; @ASB2007; @TGFGCTI2010], or to use excitations spectroscopy and photon-correlation measurements of a weakly driven system to infer certain properties of the interacting photonic system [@CGTDCI2009; @PEH2009; @KAL2010; @LH2010; @KAWLC2011].

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 [@Poyatos96; @Diehl2008; @Verstraete09; @WMLZ2010; @Kastoryano11; @Muschik11; @Barreiro11; @Krauter11]. 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 [@Diehl2008]. In the case of photonics systems, it is essential that this mechanism implements a *number-conserving* 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.

Coupled cavity arrays {#cavityarrayintro}
=====================

![\[cavity\_array\] Schematics of a coupled array of photonic cavities (represented in grey). The photons in each cavity are tunnel-coupled to neighboring sites with amplitude $J$ and decay with rate $\Gamma$. The cavities are driven with an external coherent field of strength $|\Omega|$. This external driving compensates losses and ensures a finite stationary photon population.](fig01){width="0.8\linewidth"}

Fig. \[cavity\_array\] illustrates the basic setup for an array of $L$ coupled cavities, where each cavity is represented by a single photonic mode of frequency $\omega_{\rm c}$ and bosonic annihilation operator $\hat c_\ell$. 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 Kerr-type 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. [@HBP2008]) $$\label{BHmodel}
\hat{\cal H}_{\rm c}= \sum_{\ell} \omega_{\rm c} \hat{c}_{\ell}^{\dag} \hat{c}_{\ell}  - J \sum_{\ell}  (\hat{c}_{\ell} \hat{c}_{\ell+1}^\dag +  \hat{c}_{\ell+1} \hat{c}_{\ell}^{\dag})  + \frac{U}{2}\sum_{\ell}  \hat{c}^{\dag}_{\ell} \hat{c}_{\ell}^\dag \hat{c}_{\ell} \hat{c}_{\ell}.$$ 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. (\[BHmodel\]) describes a gas of interacting bosons on a lattice, where the hopping $J$ competes with on-site interaction $U$. However, in contrast to the cold-atom physics scenario, the bosons here can decay, and under experimentally-relevant conditions $\omega_{\rm c} \gg J,U,k_{\rm B}T$ (where $T$ is the temperature and $k_{\rm 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 $\hat{\cal
  H}_{\rm 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 $\rho$ by a master equation (ME) of the form, $$\label{effectivecavity} 
\dot \rho = - i[\hat{\cal H}_{\rm c} + \hat{\cal H}_\Omega(t),\rho] + \Gamma \sum_\ell \mathcal{D}[\hat c_\ell]\rho +{\cal L}_{\kappa} \rho,$$ where $\mathcal{D}[\hat{c}]\rho \equiv 2 \hat{c} \rho \hat{c}^\dag -
\hat{c}^\dag \hat{c} \rho - \rho \hat{c}^\dag \hat{c}$. In Eq. (\[effectivecavity\]) the Hamiltonian $\hat{\cal H}_{\Omega}(t)=
\sum_{\ell} \Omega_{\ell} ( e^{-i \omega_{\rm d} t} \hat{c}_{\ell}^{\dag} + e^{i
  \omega_{\rm d} t} \hat{c}_{\ell})$ describes an external driving field of frequency $\omega_{\rm d}$ which is used to excite the system, and the second term accounts for photon losses in each cavity with a field decay rate $\Gamma$. While a finite driving field is required to counteract the losses, it will in general also compete with $\hat{\mathcal{H}}_{\rm 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 [@HBP2006; @ASB2007; @TGFGCTI2010] (where $\Omega_{\ell}=0$ for times $t>0$) [@TGFGCTI2010] or to use weak excitation spectroscopy [@CGTDCI2009; @PEH2009; @KAL2010; @LH2010; @KAWLC2011] ($\Omega_{\ell}\rightarrow 0$) to probe the many-body spectrum of Hamiltonian $\hat{\cal H}_{\rm
  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 ${\cal L}_{\kappa}$ in Eq. (\[effectivecavity\]). 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 $$\label{Lth}
{\cal L}_{\kappa} =  \sum_{\ell} \frac{\kappa}{4} \mathcal{D}[(\hat{c}_{\ell}^\dag + \hat{c}^{\dag}_{\ell+1})(\hat{c}_{\ell} - \hat{c}_{\ell+1})]  + \frac{\kappa'}{4} \mathcal{D}[(\hat{c}_{\ell}^\dag - \hat{c}_{\ell+1}^{\dag})(\hat{c}_{\ell} +\hat{c}_{\ell+1})].$$ The interpretation of this term can be seen best in the case of just two cavities. Then, the first term in Eq. (\[Lth\]) describes the scattering of photons from the asymmetric (energetically higher) mode $\hat{c}_a \equiv
(\hat{c}_1 - \hat{c}_2)/\sqrt{2}$ into the symmetric (energetically lower) mode $\hat{c}_s \equiv (\hat{c}_1 + \hat{c}_2)/\sqrt{2}$, while preserving the total photon number. The second term in Eq. (\[Lth\]) describes the reverse process. In the absence of losses, the action of ${\cal L}_{\kappa}$ 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 $\kappa'/\kappa= \exp(-2 J/k_{\rm B}T_{\rm eff})$ defining an effective temperature $T_{\rm eff}$. Similarly, for a whole cavity array, the processes $\sim \kappa$ scatter photons to lower momenta and – as previously shown for an analogous system of cold atoms [@Diehl2008] – drive the system towards a condensate of photons in the zero momentum mode, $|{\rm BEC} \rangle =
(\hat{c}_{q=0}^{\dag})^N |0\rangle / \sqrt{N!}$, where $q$ is the quasimomentum and $N$ the total number of particles. Also in this case the opposite processes $\sim \kappa'$ can be roughly interpreted as a finite temperature effect, although, as it will be shown below, for $L>2$ and $\kappa'\neq0$ the stationary state of ${\cal L}_{\kappa}$ does no longer represent a true thermal equilibrium state.

The competition of ${\cal L}_{\kappa}$ with a Kerr-type nonlinearity in the Hamiltonian has already been discussed in the atomic, number-conserving setup [@DTMFZ2010; @Tomadin11]. However, in the photonic case, under realistic conditions, ${\cal L}_{\kappa}$ competes with external driving fields and intrinsic photon losses as described by the full ME (\[effectivecavity\]). To study this competition more explicitly, we focus in the following on the limit $U \rightarrow 0$, and consider the experimental scenario where neighboring cavities are driven by coherent fields of alternating amplitude $\Omega_\ell=(-1)^\ell \Omega$ (see Fig. \[cavity\_array\]). 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 ${\cal L}_{\kappa}$.

Physical implementation
=======================

![\[circuit\] Dissipation engineering with circuit QED. a) Implementation of the Hamiltonian in Eq. (\[initialHamiltonian\]). 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. (\[rotatingframeH\]). ](fig02a "fig:"){width="0.48\linewidth"} ![\[circuit\] Dissipation engineering with circuit QED. a) Implementation of the Hamiltonian in Eq. (\[initialHamiltonian\]). 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. (\[rotatingframeH\]). ](fig02b "fig:"){width="0.48\linewidth"}

In this section we show how the photon scattering mechanism described by Eq. (\[Lth\]) can be implemented in an array of superconducting microwave resonators. The basic setting is illustrated in Fig. \[circuit\]a 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. \[circuit\]a. By restricting each resonator to a single mode $\hat{c}_{\ell}$, $\ell=1,2$, and by approximating each Cooper pair box by a two-level system with ground state $|g\rangle$ and excited state $|e\rangle$, we obtain $$\label{initialHamiltonian}
\hat{\cal H} = \sum_{\ell=1,2}  \omega_{\rm c} \hat{c}_{\ell}^\dag \hat{c}_{\ell} 
+ \sum_{\ell=1,2}  g (\hat{c}_{\ell}^\dag \hat\sigma_{-}^{(\ell)} + \hat{c}_{\ell} \hat\sigma_{+}^{(\ell)}) +  \hat{\cal H}_{\rm q} (t),$$ where $\hat\sigma_{-}^{(\ell)} \equiv | g_{\ell} \rangle \langle e_{\ell} |$. The first part is the free cavity Hamiltonian, while the second term describes the dipole interaction between photons and qubits with strength $g$. Finally, $\hat{\cal H}_{\rm q}(t)$ is the Hamiltonian of the two coupled qubits, which we assume of the generic form [^1] $$\begin{aligned}
 \label{qubitH}
\hat{\cal H}_{\rm q}(t)&=& \sum_{\ell=1,2}[ \omega_{\rm q} |e_{\ell} \rangle \langle e_{\ell}| + \Omega_{\rm q}^{(\ell)} ( e^{-i\omega_0 t} \hat\sigma_+^{(\ell)} + e^{i\omega_0 t} \hat\sigma_{-}^{(\ell)}) ]\nonumber \\
& &+ \lambda_z \hat\sigma_{ee}^{(1)} \hat\sigma_{ee}^{(2)} + \lambda_x (\hat\sigma_{+}^{(1)} \hat\sigma_{-}^{(2)} + \hat\sigma_{-}^{(1)} \hat\sigma_{+}^{(2)}).\end{aligned}$$ Here, $\Omega_{\rm q}^{(\ell)}$ are the Rabi frequencies of an applied microwave field, which is used to drive the qubits at frequency $\omega_0$. The coupling constants $\lambda_z$ and $\lambda_x$ are the strengths of the diagonal and off-diagonal qubit-qubit interactions, respectively. The Hamiltonian (\[initialHamiltonian\]) can be rewritten in terms of the symmetric $\hat c_s \equiv (\hat{c}_1 +
\hat{c}_2)/\sqrt{2}$ and antisymmetric $\hat c_a\equiv (\hat{c}_1 - \hat{c}_2)/\sqrt{2}$ cavity modes, and the qubit states $|E\rangle\equiv |ee\rangle$, $|G\rangle\equiv |gg\rangle$, $|S\rangle\equiv (|eg\rangle+ |ge\rangle)/\sqrt{2}$, and $|A\rangle\equiv
(|eg\rangle- |ge\rangle)/\sqrt{2}$. In this new basis, choosing $\Omega_{\rm
  q}^{(1)} = -\Omega_{\rm q}^{(2)} \equiv \Omega_a/\sqrt{2}$, and changing into a frame rotating with $\omega_0$, we obtain $$\begin{aligned}
 \label{rotatingframeH}
\hat{\cal H} & = & -\delta \hat{c}_s^\dag \hat{c}_s - \delta \hat{c}_a^\dag \hat{c}_a  - \Delta_e  |E\rangle\langle E|  - \Delta_s  |S\rangle\langle S| - \Delta_a |A\rangle\langle A| \nonumber \\
&& + \Omega_a [ |A\rangle \langle G| -  |E\rangle \langle A| + {\rm H.c.}] \nonumber \\
& & + g [ \hat{c}_s^\dag (  |G\rangle \langle S| + |S\rangle \langle E| ) + \hat{c}_a^\dag ( |G\rangle \langle A| - |A\rangle \langle E| ) + {\rm H.c.}],\end{aligned}$$ where $\delta= \omega_0-\omega_{\rm c}$, $\Delta_{s,a} = \omega_0 - \omega_{\rm
  q} \mp \lambda_x$, $\Delta_e = 2(\omega_0 - \omega_{\rm q})-\lambda_z$. Apart from the coherent dynamics described by Eq. (\[rotatingframeH\]), we include dissipation in the form of intrinsic cavity loss with rate $2\Gamma$ and qubit decay with rate $\gamma$. The full system dynamics is then described by the ME $$\label{CircuitQEDME}
\dot \rho = -i[\hat{\cal H},\rho] + \sum_{\eta = s,a} \left (   \Gamma  \mathcal{D}[\hat{c}_\eta] \rho + \frac{\gamma}{2} \mathcal{D}[|G\rangle\langle \eta|] \rho + \frac{\gamma}{2} \mathcal{D}[|\eta \rangle\langle E|] \rho \right ).$$ and a summary of the energy levels and the most relevant transitions is shown in Fig. \[circuit\]b.

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\rangle$ to $|S\rangle$ by absorbing a photon from the antisymmetric cavity mode ($\hat{c}_a$) and emitting a photon into the symmetric cavity mode ($\hat{c}_s^{\dagger}$) (see Fig. \[circuit\]b). Since the overall process is $\mathcal{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. \[circuit\]b. In particular, we assume that $\delta$, $\Delta_e$, and $\Delta \equiv \delta-\Delta_e$ are much larger than all the other frequency scales $g, \Omega_a, \Delta_{s,a},
\Gamma, \gamma$. 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 transformation $\tilde{\cal H}
\equiv \hat{V} \hat{\cal H} \hat{V}^\dag$ with $\hat{V}=e^{\hat{S}}$, and make the ansatz $$\hat{S} = \hat{c}_s^\dag  \left [ \vphantom{\sum}  \alpha_{1,s} |G\rangle \langle S| + \alpha_{2,s} |S\rangle \langle E| \right ]  + \hat{c}_a^\dag \left [ \vphantom{\sum} \alpha_{1,a} |G\rangle \langle A| - \alpha_{2,a} |A\rangle \langle E| \right ] - {\rm H.c.}~$$ We define $\hat{\cal H}_0$ and $\hat{\cal H}_g$ as the first and last lines of the Hamiltonian (\[rotatingframeH\]), respectively, and choose the coefficients $\alpha_i$ such that $[\hat{S},\hat{\cal H}_0] = -\hat{\cal H}_g$. This can be achieved by setting $$\alpha_{1,\eta} = \frac{g}{\delta-\Delta_\eta},
\qquad \alpha_{2,\eta}= \frac{g}{\delta-(\Delta_e -\Delta_\eta)},
\qquad \eta = s, a.$$ In view of the assumptions discussed above, we will use the approximate results $\alpha_{1,\eta}\approx g/\delta$, $\alpha_{2,\eta}\approx g/\Delta$ and define the parameter $\epsilon=\Delta/\delta <1$. After this transformation the Hamiltonian (\[rotatingframeH\]) reads $$\tilde{\cal H} =  \tilde{\cal H}_{\rm q} + \tilde{\cal H}_{\rm c} + \tilde{\cal H}_{\rm int} + \mathcal{O}(g^3, g \Omega_a).$$ Here $\tilde{\cal H}_{\rm q}$ is the qubit Hamiltonian $\hat{\cal H}_{\rm q}$ written above, with small frequency shifts absorbed into a redefinition of the detunings $\Delta_{e,s,a}$. The modified cavity Hamiltonian is $$\tilde{\cal H}_{\rm c} = - \bar \delta \hat{c}_s^{\dag} \hat{c}_s - \bar \delta \hat{c}_a^{\dag} \hat{c}_a - J (\hat{c}_s^{\dag} \hat{c}_s - \hat{c}_a^{\dag} \hat{c}_a) + g_{\rm eff} \left[( \hat{c}_s^\dag \hat{c}_a+ \epsilon \hat{c}_a^\dag \hat{c}_s )  \langle P_{sa} \rangle + {\rm H.c.} \right],$$ where $g_{\rm eff}\equiv g^2/\Delta$, $\bar \delta \equiv \delta - g_{\rm eff}
(\langle P_{aa} \rangle + \langle P_{ss} \rangle)/2$, $J\equiv g_{\rm eff}
(\langle P_{aa} \rangle- \langle P_{ss} \rangle)/2$, $P_{sa} \equiv
|S\rangle\langle A|$, $$\begin{aligned}
P_{ss} & = &   (|S\rangle\langle S| - |E\rangle\langle E|) + \epsilon  (|G\rangle\langle G| -|S\rangle\langle S|), \\
P_{aa} & = &  (|A\rangle\langle A| - |E\rangle\langle E|) + \epsilon  (|G\rangle\langle G| -|A\rangle\langle A|),\end{aligned}$$ and the average $\langle \cdot \rangle$ is taken with respect to the stationary qubit state in the limit $g\rightarrow 0$. For the parameter regime of interest, $\langle P_{sa} \rangle\approx 0$, and $\tilde{\cal H}_{\rm 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 $$\tilde{\cal H}_{\rm int} \simeq  g_{\rm eff} \left [ \hat{c}_s^\dag \hat{c}_s \bar P_{ss} + \hat{c}_a^\dag \hat{c}_a  \bar P_{aa} +(\hat{c}_s^\dag \hat{c}_a+ \epsilon \hat{c}_a^\dag \hat{c}_s)  \bar P_{sa} + (\hat{c}_s^\dag \hat{c}_a+ \epsilon \hat{c}_a^\dag \hat{c}_s)^\dag \bar P_{sa}^\dag \right ],~$$ where $\bar P_{km}=P_{km} -\langle P_{km}\rangle$. Note that in $\tilde{\cal H}_{\rm
  int}$ we have only kept resonant two-photon processes and already omitted terms like $\sim g\Omega_a/\Delta \times \hat{c}_{s,a}$ and $g^2/\Delta \times
\hat{c}_{s,a}^2$. Both of these terms oscillate with the detuning $\delta$ and will only give small corrections to the results presented below.

The transformation $\hat{V}$ induces a weak mixing between photon and qubit states, which apart from generating new effective interactions also modifies the dissipative couplings. Since we consider $\gamma \gg \Gamma $, we find that the main result of this mixing is an enhancement of the cavity decay rates $ \Gamma
\rightarrow \Gamma_{a,s} \equiv \Gamma + \tilde\Gamma_{s,a}$. After a transformation of the jump operators in the ME (\[CircuitQEDME\]), and averaging over the qubit states, we find that for the parameter regime of interest, these rates are $$\begin{aligned}
\tilde\Gamma_{s}&\approx& \frac{\gamma g^2}{2\Delta^2} \epsilon^2 \left[ \langle |G\rangle\langle G| \rangle + \langle | A \rangle \langle A | \rangle \right], \\
\tilde\Gamma_{a}&\approx& \frac{\gamma g^2}{2\Delta^2} \left[ \epsilon^2 \langle |G\rangle\langle G| \rangle+ (2+\epsilon^2) \langle |A\rangle\langle A| \rangle \right].\end{aligned}$$

In summary, we find that the system dynamics in the new dressed state representation is well described by the ME $$\dot \rho = ( \mathcal{L}_{\rm q} +  \mathcal{L}_{\rm c}  + \mathcal{L}_{\rm int} ) \rho,$$ where $$\begin{aligned}
\mathcal{L}_{\rm q} \rho = -i[\tilde{\cal H}_{\rm q},\rho] + \sum_{\eta=s,a} \left ( \frac{\gamma}{2} \mathcal{D}[|G\rangle\langle \eta|] \rho + \frac{\gamma}{2} \mathcal{D}[|\eta \rangle\langle E|] \rho \right ),\\
\mathcal{L}_{\rm c} \rho = -i [\tilde{\cal H}_{\rm c},\rho] + \sum_{\eta=s,a} \Gamma_{\eta} \mathcal{D}[\hat{c}_{\eta}] \rho,\\
\mathcal{L}_{\rm int}\rho = -i[\tilde{\cal H}_{\rm int},\rho].\end{aligned}$$ In the limit $g_{\rm eff}\rightarrow 0 $, qubits and cavities are decoupled, and the system relaxes on a timescale $\gamma^{-1}$ into the state $\rho(t) \simeq \rho_{\rm c} (t) \otimes \rho_{\rm q}^0$, where $\mathcal{L}_{\rm q} \rho_{\rm q}^0 = 0$. For finite $g_{\rm eff} <\gamma$, 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 [@BP2002], $$\dot \rho_{\rm c} = {\cal L}_{\rm c} \rho - \int_0^\infty d\tau \, {\rm Tr}_{\rm q}\{ [\tilde{\cal H}_{\rm int}, e^{\tau \mathcal{L}_{\rm q}} \left( [\tilde{\cal H}_{\rm int}, \rho_{\rm c}(t) \otimes \rho_{\rm q}^0] \right) ] \}.$$ Evaluating this expression we obtain $$\begin{aligned}
\label{eq:EffCavityME}
\dot \rho_{\rm c} & \approx & -i[\tilde{\cal H}_{\rm c} + \hat{\cal H}_{\rm eff},\rho_{\rm c} ] + \sum_{\eta=s,a} \Gamma_{\eta} \mathcal{D}[\hat{c}_{\eta}] \rho_{\rm c} +   \kappa \mathcal{D}[(\hat{c}_s^\dag \hat{c}_a+ \epsilon \hat{c}_a^\dag \hat{c}_s)] \rho_{\rm c}  \nonumber \\
& &+ \sum_{\eta,\eta'=s,a}\Gamma^\Phi_{\eta\eta'} \left(\hat  n_{\eta} \rho_{\rm c} \hat n_{\eta'} + \hat n_{\eta'} \rho_{\rm c} \hat n_{\eta}
 - \hat n_{\eta} \hat n_{\eta'} \rho_{\rm c} - \rho_{\rm c} \hat n_{\eta} \hat n_{\eta'} \right ),\end{aligned}$$ where $\hat n_\eta=\hat c_\eta^\dag \hat c_\eta$ and $\hat{\cal H}_{\rm eff}$ is an effective photon-photon interaction $$\hat{\cal H}_{\rm eff}=  U_{\rm eff} (\hat{c}_s \hat{c}_a^\dag +\epsilon \hat{c}_a \hat{c}^\dag_s)(\hat{c}_s^\dag \hat{c}_a +\epsilon \hat{c}_a^\dag \hat{c}_s) +  \sum_{\eta,\eta'=s,a} U_{\eta\eta'} \hat n_{\eta} \hat n_{\eta'}.$$ In Eq. (\[eq:EffCavityME\]) we have defined $\kappa \equiv g_{\rm eff}^2 {\rm
  Re} \{S_{as,sa}(0)\}$, $\Gamma^\Phi_{\eta\eta'} \equiv g_{\rm eff}^2 {\rm Re}
\{S_{\eta\eta,\eta'\eta'}(0)\}$, and the interactions $U_{\rm eff} \equiv
g_{\rm eff}^2 {\rm Im} \{S_{as,sa}(0)\}$ and $U^\Phi_{\eta\eta'} \equiv g_{\rm
  eff}^2 {\rm Im} \{S_{\eta\eta,\eta'\eta'}(0) \} $, given in terms of the qubit correlation spectra $$\label{spectrumEq}
S_{\eta\eta',\varsigma\varsigma'}(\omega) \equiv  \int_0^\infty d\tau \, e^{-i\omega \tau} \langle \bar P_{\eta\eta'}(\tau) \bar P_{\varsigma\varsigma'}(0)\rangle,$$ which can be calculated using the quantum regression theorem [@GardinerZollerBook].

Discussion
----------

By ignoring coherent photon-photon interactions for the moment, we see that the first line of the effective cavity dynamics given in Eq. (\[eq:EffCavityME\]) represents the dissipatively coupled cavities as introduced in Eq. (\[effectivecavity\]). In particular, we can write $$\label{expandLth}
\kappa \mathcal{D}[(\hat{c}_s^\dag \hat{c}_a+ \epsilon \hat{c}_a^\dag \hat{c}_s)]  = \kappa  \mathcal{D}[\hat{c}_s^\dag \hat{c}_a] + \kappa \epsilon^2  \mathcal{D}[\hat{c}_a^\dag \hat{c}_s]+ \mathcal{L}_\epsilon.$$ By comparison of Eqs. (\[eq:EffCavityME\]) and (\[Lth\]) we identify $\kappa'=\epsilon^2
\kappa$, with the parameter $\epsilon^2\equiv (\Delta/\delta)^2$ determining the minimal effective temperature of the engineered two-photon process. For typical parameters, $\omega_{\rm q}=10$ GHz, $\delta=1$ GHz and $\Delta_e=0.8$ GHz, we obtain $\epsilon=0.2$, and hence the limit $\kappa'\rightarrow 0$ considered in most parts of the paper can be implemented to a very good approximation. However, we point out that the ratio $\kappa'/\kappa$ can always be increased artificially, for example, by adding an additional coherent or incoherent driving field to populate the state $|S\rangle$. Our derivation also shows the existence of coherences between $\hat{c}_s^\dag \hat{c}_a$ and $\hat{c}_a^\dag \hat{c}_s$, which lead to additional squeezing terms of order $\mathcal{O}(\kappa \epsilon)$, and that in Eq. (\[expandLth\]) are summarized by $ \mathcal{L}_\epsilon$. For $\epsilon \ll1 $, 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 $\kappa'$. Indeed, numerical simulations similar to those discussed in Sec. \[Results\] show no significant differences between the exact ME (\[eq:EffCavityME\]) and our model Liouvillian (\[Lth\]), even up to values of $\epsilon\approx 0.5$. However, by choosing $\epsilon\approx 1$ and $\Gamma\ll \kappa$, it would be possible to extend our model and study physical effects directly linked to the parameter $\epsilon$, such as number-conserving squeezing as previously analyzed in [@WM2011].

![\[FOM\] a) Single and two-photon rates which appear in the effective model of Eq. (\[eq:EffCavityME\]). The photon scattering rate of interest, $\kappa$, becomes significantly larger than the others in a certain region of the parameter space. Here we have chosen $g=25$ MHz, $\delta=1$ GHz, $\Delta_s=\Delta_a=0$ MHz, $\Delta_e=800$ MHz, $\gamma=10$ MHz, $
  \Gamma =10$ kHz. b) Figure of merit for the validity of our model for high-momentum modes, $F_a\equiv\kappa/(\Gamma_a+\Gamma_{aa}^\Phi)$, as a function of qubit driving $\Omega_a$ and detuning $\Delta_a$. Rest of parameters as in (a). c) The same plot as in (b) but for low-momentum modes, $F_s\equiv\kappa/(\Gamma_s+\Gamma_{ss}^\Phi)$. ](fig03a "fig:"){width="0.48\linewidth"} ![\[FOM\] a) Single and two-photon rates which appear in the effective model of Eq. (\[eq:EffCavityME\]). The photon scattering rate of interest, $\kappa$, becomes significantly larger than the others in a certain region of the parameter space. Here we have chosen $g=25$ MHz, $\delta=1$ GHz, $\Delta_s=\Delta_a=0$ MHz, $\Delta_e=800$ MHz, $\gamma=10$ MHz, $
  \Gamma =10$ kHz. b) Figure of merit for the validity of our model for high-momentum modes, $F_a\equiv\kappa/(\Gamma_a+\Gamma_{aa}^\Phi)$, as a function of qubit driving $\Omega_a$ and detuning $\Delta_a$. Rest of parameters as in (a). c) The same plot as in (b) but for low-momentum modes, $F_s\equiv\kappa/(\Gamma_s+\Gamma_{ss}^\Phi)$. ](fig03bc "fig:"){width="0.48\linewidth"}

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 $\Gamma_{s,a}$ introduced above, and on the other hand, we obtain additional dephasing terms $\Gamma_{ij}^\Phi$ due to fluctuating Stark shifts when the qubit jumps incoherently between states $|G\rangle$ and $|A\rangle$. Since the ratio $\kappa/\Gamma_{s,a}$ scales as $\sim
g^2/\gamma^2$, single photon losses can be suppressed in the strong coupling regime, as long as $\kappa$ also exceeds the bare decay rate $ \Gamma $. The ratio $\kappa/\Gamma_{ij}^\Phi$ is independent of $g$, and depends on the details of the correlation functions (\[spectrumEq\]) and therefore on the parameters $\Omega_a$, $\Delta_{s,a}$. In Fig. \[FOM\]a, $\kappa$ is compared with the other decoherence rates, and we see that an experimentally feasible parameter range exists, where $\kappa$ is the dominant rate. The main corrections arise from the enhanced decay $\Gamma_a$ and the dephasing $\Gamma_{aa}^\Phi$ of the antisymmetric mode, while the corresponding rates for the symmetric mode are much smaller. In Fig. \[FOM\]b and \[FOM\]c we plot the two quantities $F_{\eta}\equiv \kappa / (\Gamma_\eta +
\Gamma_{\eta\eta}^\Phi)$, $\eta=s,a$, as a function of the control parameters $\Omega_a$, $\Delta_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 $\kappa$ can exceed all decoherence processes discussed here, and therefore the proposed model in Eq. (\[effectivecavity\]) 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 $\hat{\cal H}_{\rm eff}$. For a resonant two-photon process, $\Delta_s=0$, we find $U_{\rm eff}\approx 0$. Moreover, $U_{\eta\eta'}^\Phi\approx \Gamma_{\eta\eta'}^\Phi < \kappa$ for the typical parameters considered above, and in particular these interactions vanish for $\Delta_a=0$. Therefore, in this work we focus exclusively on the purely dissipative cavity dynamics. However, by setting $|\Delta_s|>\gamma$ instead, the regime of strong coherent photon-photon coupling $U_{\rm
  eff}>\kappa>\Gamma_{\eta\eta'}^\Phi,\Gamma_{s,a}$ can be engineered as well, having $\hat{\cal H}_{\rm eff}\simeq U_{\rm eff} \hat{c}_s^\dag \hat{c}_s
\hat{c}_a^\dag \hat{c}_a$.

In summary, we find that the ME (\[effectivecavity\]) 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 $\lambda_{x,z}$ plays the role of the nonlinearity.

Condensation of photons: two cavities {#Results}
=====================================

Semiclassical approximation
---------------------------

In the previous section we have shown how to implement the effective ME (\[effectivecavity\]). Now, we discuss how the scattering between photons induced by the Liouvillian ${\cal L}_{\kappa}$ 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. $\Omega_1=-\Omega_2=\Omega/\sqrt{2}$. In the rotating frame with respect to the drive frequency $\omega_{\rm d}$ the equations of motion (EOM) for the field expectation values $\langle\hat{c}_{s}\rangle$ and $\langle\hat{c}_{a}\rangle$ of the symmetric and antisymmetric cavity modes are $$\begin{array}{l} \label{EOMcoherentpumping}
\partial_{t} \langle \hat{c}_s \rangle = [- i (\delta_{\rm d}-2J) - \Gamma] \langle \hat{c}_s \rangle + \kappa \langle \hat{c}_s \hat{c}_a^{\dag} \hat{c}_a \rangle,\\
\partial_{t} \langle \hat{c}_a \rangle = [- i \delta_{\rm d} - \Gamma ] \langle \hat{c}_a \rangle - i \Omega - \kappa \langle (\hat{c}_s^{\dag} \hat{c}_s + 1) \hat{c}_a \rangle,
\end{array}$$ where $\delta_{\rm d}=\omega_{\rm c}+J-\omega_{\rm d}$ is the detuning of the antisymmetric mode with respect to the driving field. For the moment, we assume $\kappa' = 0$ in Eq. (\[Lth\]) 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 $$\begin{array}{l} \label{meanfieldEOM2cavities}
\partial_{t} s = - i (\delta_{\rm d}-2J) s - (\Gamma - \kappa |a|^{2})s, \\
\partial_{t} a = - i \delta_{\rm d} a - i  \Omega - [\Gamma + \kappa(|s|^{2} + 1)] a,
\end{array}$$ where $s \equiv \langle \hat{c}_s \rangle$ and $a \equiv \langle \hat{c}_a
\rangle$ 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 (\[EOMcoherentpumping\]).

![\[transitionmeanfield\] Results (\[semiclassicalresults\]) of the semiclassical approximation to the EOM (\[EOMcoherentpumping\]). a) Population fraction $n_s/n$ in the symmetric mode as a function of cavity pumping $\Omega$ and dissipation strength $\kappa$. b) Population in the symmetric (solid line) and antisymmetric (dashed line) mode for $\kappa=\Gamma$ as a function of the pumping strength $\Omega$. The threshold (\[threshold\]) is $\Omega_{\rm crit} = 2.0 \, \Gamma$. ](fig04a "fig:"){width="0.48\linewidth"} ![\[transitionmeanfield\] Results (\[semiclassicalresults\]) of the semiclassical approximation to the EOM (\[EOMcoherentpumping\]). a) Population fraction $n_s/n$ in the symmetric mode as a function of cavity pumping $\Omega$ and dissipation strength $\kappa$. b) Population in the symmetric (solid line) and antisymmetric (dashed line) mode for $\kappa=\Gamma$ as a function of the pumping strength $\Omega$. The threshold (\[threshold\]) is $\Omega_{\rm crit} = 2.0 \, \Gamma$. ](fig04b "fig:"){width="0.48\linewidth"}

The exact equations (\[EOMcoherentpumping\]) 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 [@HakenBook], 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 $$\label{threshold}
\Omega_{\rm crit} \equiv \Gamma\sqrt{\frac{\delta_{\rm d}^2 + (\Gamma+\kappa)^2}{\Gamma\kappa}},$$ for the pumping strength $\Omega$. Below the threshold, the only stable steady-state solution to Eq. (\[meanfieldEOM2cavities\]) is (c.f. also \[AppendixB\]) $$\begin{aligned}
 \label{solbelowthreshold}
|s|^2=0, \qquad |a|^2= \frac{\Omega^2}{\delta_{\rm d}^2+(\Gamma+\kappa)^2},\end{aligned}$$ and only the pumped antisymmetric mode is populated. For driving strengths above the threshold we obtain the stable solution $$\begin{aligned}
 \label{solabovethreshold}
|s|^2=\sqrt{\frac{\Omega^2}{\Gamma\kappa} -\frac{\delta_{\rm d}^2}{\kappa^2}} - \frac{\Gamma}{\kappa} - 1, \;\;\;\;\; |a|^2= \frac{\Omega^2}{\delta_{\rm d}^2+[\Gamma+\kappa(|s|^2+1)]^2},\end{aligned}$$ where the symmetric mode features nonvanishing occupation. Note that these equations depend on $J$ only via the detuning $\delta_{\rm d}$, which has the only effect to reduce the driving strength. Therefore, from now on we will restrict ourselves to $\delta_{\rm d}=0$ and $J\rightarrow 0$. The stationary solution to the semiclassical EOM then simplifies to $$\begin{array}{l} \label{semiclassicalresults}
|s|^2 = 0, \qquad |a|^2=\frac{\Omega^2}{(\Gamma+\kappa)^2}, \qquad \mbox{if} \quad \Omega < \Omega_{\rm crit}, \\
|s|^2 = \frac{\Omega}{\sqrt{\Gamma\kappa}} - \frac{\Gamma}{\kappa} - 1, \qquad |a|^2=\frac{\Gamma}{\kappa}, \qquad \mbox{if} \quad \Omega > \Omega_{\rm crit}.
\end{array}$$ The behavior of the populations $n_s\equiv |s|^2$ and $n_a\equiv |a|^2$ is plotted in Fig. \[transitionmeanfield\] as a function of the pumping intensity and the effective decay $\kappa$. Below threshold, the symmetric mode is empty and the antisymmetric mode occupation grows quadratically with the pumping intensity. For $\Omega>\Omega_{\rm crit}$, the antisymmetric mode population saturates, and the symmetric mode occupation grows linearly. Thus, if $\Omega$ 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 (\[Lth\]).

The discussion of Eqs. (\[semiclassicalresults\]) illustrates both, similarities and differences between the present setup and a laser [@AINYPNT2007; @RIA2007; @YLSN2007; @HFHSS2008]. 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 ${\cal L}_{\kappa}$ 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 [@Snoke2003; @SL2002; @SLS2003; @Scully1999; @KSZZ2000] 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 $\kappa\to 0$ as well as for $\kappa \rightarrow \infty$ the condensation is suppressed, and that the phase diagram in Fig. \[transitionmeanfield\] is symmetric with respect to the ratio $\kappa/\Gamma$. For small $\kappa$, not enough photons are scattered into the symmetric mode, while for large $\kappa$ the antisymmetric mode is highly overdamped and the population of the whole system is small as in the first place.

Exact diagonalization
---------------------

![\[numerical2cavities\] Numerical solution to the exact equations of motion for $L=2$ cavities. The occupation numbers $n_s$ and $n_a$ of the symmetric and antisymmetric mode are plotted as a function of the coherent driving strength $\Omega$ for a) $\kappa/\Gamma=0.2$ and b) $\kappa/\Gamma=5$. In both plots $\delta_{\rm d}=0$, and the dashed lines indicate the corresponding semiclassical predictions given by Eq. (\[semiclassicalresults\]). The dotted lines show the two-photon correlation function of the symmetric mode $g^{(2)}=\langle
\hat{n}_s (\hat{n}_s - 1) \rangle / \langle \hat{n}_s \rangle^{2}$. ](fig05a "fig:"){width="0.48\linewidth"} ![\[numerical2cavities\] Numerical solution to the exact equations of motion for $L=2$ cavities. The occupation numbers $n_s$ and $n_a$ of the symmetric and antisymmetric mode are plotted as a function of the coherent driving strength $\Omega$ for a) $\kappa/\Gamma=0.2$ and b) $\kappa/\Gamma=5$. In both plots $\delta_{\rm d}=0$, and the dashed lines indicate the corresponding semiclassical predictions given by Eq. (\[semiclassicalresults\]). The dotted lines show the two-photon correlation function of the symmetric mode $g^{(2)}=\langle
\hat{n}_s (\hat{n}_s - 1) \rangle / \langle \hat{n}_s \rangle^{2}$. ](fig05b "fig:"){width="0.48\linewidth"}

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}=\langle \hat c_s^\dag \hat c_s\rangle $ and $n_{a}=\langle \hat c_a^\dag \hat c_a\rangle$ in the symmetric and antisymmetric mode, respectively, are shown in Fig. \[numerical2cavities\] as a function of the pumping strength $\Omega$ and for two different values of $\kappa/\Gamma$. In the limit $\kappa\ll\Gamma$, 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 $\Omega_{\rm crit}$ the population $n_{s}$ grows linearly with the driving strength, while the antisymmetric mode saturates at a value $n_a=\Gamma/\kappa$. 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 two-photon correlation $g^{(2)} \equiv \langle
\hat{n}_s (\hat{n}_s - 1) \rangle / \langle \hat{n}_s \rangle^{2}$ reveals that, when crossing $\Omega_{\rm crit}$, the photon statistics of the symmetric mode changes from a thermal ($g^{(2)}\approx 2$) to a Poissonian ($g^{(2)}\approx 1$) distribution. This aspect is in agreement with the standard lasing transition.

In the opposite regime, $\kappa\gg\Gamma$, Eq. (\[semiclassicalresults\]) 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. \[numerical2cavities\]b shows that in this limit the transition is completely washed out and $n_s$ exceeds $n_a$ for all pumping strengths. Nevertheless, above $\Omega_{\rm crit}$ we still observe the characteristic saturation of $n_a$ and – apart from an additional offset – the 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 $\hat c_a$ mode is given by $\Gamma +\kappa(\langle \hat n_s\rangle +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’ $\hat c_a$. The limit $\kappa\gg \Gamma$ and $\Omega\approx \Omega_{\rm crit}$ would then result in a ‘single photon condensate’ with $ n_a\ll n_s \sim 1$ and a non-classical statistics with $g^{(2)}\rightarrow 0.5$.

Effective temperature
---------------------

Finally, let us briefly comment on the effect of the reverse photon scattering process $\sim \kappa'$ in Eq. (\[Lth\]). As we have discussed in Sec. \[cavityarrayintro\], 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 $\kappa'=\kappa e^{-2 J/k_BT_{\rm
    eff}}$. By including this term, the EOM are $$\begin{array}{l}
\partial_{t} \langle \hat{c}_s \rangle = - (\Gamma+\kappa') \langle \hat{c}_s \rangle + (\kappa-\kappa') \langle \hat{c}_s \hat{c}_a^{\dag} \hat{c}_a \rangle, \\
\partial_{t} \langle \hat{c}_a \rangle = - (\Gamma+\kappa') \langle \hat{c}_a \rangle - i  \Omega - (\kappa-\kappa') \langle (\hat{c}_s^{\dag} \hat{c}_s + 1) \hat{c}_a \rangle.
\end{array}$$ These equations are of the same form of (\[EOMcoherentpumping\]), with the replacements $\Gamma\to\Gamma+\kappa'$ and $\kappa\to\kappa-\kappa'$. Therefore, in the semiclassical regime we obtain the same physics, but with a renormalized pumping threshold $$\label{finiteTthreshold}
\Omega_{\rm crit} (T_{\rm eff})\equiv (\Gamma+\kappa) \sqrt{\frac{(\Gamma/\kappa + e^{-2 J/k_BT_{\rm eff}}  )}{(1-e^{-2 J/k_BT_{\rm eff}} )}} .$$ In particular, for any fixed driving strength we can use Eq. (\[finiteTthreshold\]) to define a critical effective temperature above which photon condensation does not occur.

Photon condensation in a cavity array {#cavityarraysec}
=====================================

We turn now to the solution of Eq. (\[effectivecavity\]) in a lattice with $L
\gg 1$ cavities. The central goal of this section is to show that the many-body Liouvillian (\[Lth\]) 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 $\hat{c}_{\ell}$ in the array, with $\ell = 1, \dots L$, is coherently driven by a classical field with staggered amplitude $\Omega_{\ell} = (-1)^{\ell} \Omega$. 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. (\[effectivecavity\]) in terms of the annihilation operators $\hat{c}_q = \frac{1}{\sqrt{L}} \sum_\ell e^{i q \ell} \hat{c}_\ell$ of the photonic modes in momentum space, where $q$ takes values $-\pi (L-2)/L,\dots,\pi
(L-2)/L,\pi$ (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 by $$\label{Hamiltonianq}
\hat{\cal H}_{\rm c} + \hat{\cal H}_{\Omega} = \sum_{q} [\delta_{\rm d}-2J(1+\cos(q))] \hat{c}_{q}^{\dag} \hat{c}_{q} + \sum_{q} ( \Omega_{q} \hat{c}_{q}^{\dag} + \Omega_{q}^{\ast} \hat{c}_{q}),$$ where $\Omega_q
= \frac{1}{\sqrt{L}} \sum_{\ell} e^{iq\ell}\Omega_{\ell}$ and $\delta_{\rm d}= \omega_{\rm c}+2J-\omega_{\rm d}$ is the detuning of the $q=\pi$ mode from the driving frequency. The Liouville operators are $$\label{LiouvillianGammaq}
\Gamma \sum_{\ell} {\cal D}[\hat{c}_{\ell}] \rho = \Gamma \sum_{q} [2 \hat{c}_{q} \rho \hat{c}_{q}^{\dag} - \hat{c}_{q}^{\dag} \hat{c}_{q} \rho - \rho \hat{c}_{q}^{\dag} \hat{c}_{q}],$$ and $$\begin{aligned}
 \label{Liouvilliankappaq}
{\cal L}_{\kappa} \rho & = & \sum_{q_{1}q_{2}q_{3}q_{4}} \delta(q_{1}+q_{3}-q_{2}-q_{4}) K(q_{1},q_{2},q_{3},q_{4}) \nonumber \\
& & \times [2 \hat{c}_{q_{1}}^{\dag} \hat{c}_{q_{2}} \rho \hat{c}_{q_{3}}^{\dag} \hat{c}_{q_{4}} - \hat{c}_{q_{3}}^{\dag} \hat{c}_{q_{4}} \hat{c}_{q_{1}}^{\dag} \hat{c}_{q_{2}} \rho - \rho \hat{c}_{q_{3}}^{\dag} \hat{c}_{q_{4}} \hat{c}_{q_{1}}^{\dag} \hat{c}_{q_{2}}].\end{aligned}$$ In the latter equation, the function $K$ plays the role of a scattering kernel, which has the general form $$\begin{aligned}
K(q_{1},q_{2},q_{3},q_{4}) =  \frac{4}{L} \Big[ &&\kappa \cos\frac{q_{1}}{2} \sin\frac{q_{2}}{2} \sin\frac{q_{3}}{2} \cos\frac{q_{4}}{2} \nonumber \\
&&+\kappa'\sin\frac{q_{1}}{2} \cos\frac{q_{2}}{2} \cos\frac{q_{3}}{2} \sin\frac{q_{4}}{2}\Big].\end{aligned}$$ In the following analysis, we will focus on the case $\kappa'=0$. Furthermore, as discussed above, the photon condensation effect does not depend on the tunneling amplitude and we can assume $J \rightarrow 0$ in the Hamiltonian.

Solving the ME (\[effectivecavity\]) 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. \[condensation\_fintite\_q\]b 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. \[Results\] 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
---------------------------

![\[condensation\_fintite\_q\] a) Numerical solution of the semiclassical equations of motion (\[eq:boltzmanncoherent\]). $L=22$ momentum states are used in the integration and the population $n_{q} =
  |\psi_{q}|^{2}$ is shown in the positive half $q>0$ of the Brillouin zone. b) One-particle density matrix $\langle \hat{c}^{\dag}_{q} \hat{c}_{q'} \rangle$ in momentum space, obtained with the numerical integration of the exact equations of motion for $L=4$ cavities, integrated until $t_{\rm max} \,
  \Gamma = 0.5$. In both panels we use $\kappa = 5.0 \, \Gamma$ and $\Omega_{q=\pi} = 3.5 \sqrt{L} \, \Gamma$. ](fig06a "fig:"){width="0.48\linewidth"} ![\[condensation\_fintite\_q\] a) Numerical solution of the semiclassical equations of motion (\[eq:boltzmanncoherent\]). $L=22$ momentum states are used in the integration and the population $n_{q} =
  |\psi_{q}|^{2}$ is shown in the positive half $q>0$ of the Brillouin zone. b) One-particle density matrix $\langle \hat{c}^{\dag}_{q} \hat{c}_{q'} \rangle$ in momentum space, obtained with the numerical integration of the exact equations of motion for $L=4$ cavities, integrated until $t_{\rm max} \,
  \Gamma = 0.5$. In both panels we use $\kappa = 5.0 \, \Gamma$ and $\Omega_{q=\pi} = 3.5 \sqrt{L} \, \Gamma$. ](fig06b "fig:"){width="0.48\linewidth"}

For many cavities, the EOM for the cavity fields $\psi_{q} = \langle \hat{c}_{q} \rangle$ (analogous to Eq. (\[meanfieldEOM2cavities\])) read $$\begin{aligned}
\label{eq:boltzmanncoherent}
\partial_{t} \psi_{q} & = & [- i \delta_{\rm d} - \Gamma - \kappa (1 - \cos{q})] \psi_{q} - i \Omega_{q} \nonumber \\
& & + \frac{ \kappa}{L} \sum_{q_{1}q_{2}} \psi_{q_{1}} \psi_{q_{2}} \psi_{q_{1} + q_{2} - q}^{\ast}
[\cos{q} - \cos{(q_{1} + q_{2} - q)}],\end{aligned}$$ where $\Omega_{q} = \delta_{q,\pi} \Omega_\pi $. The numerical solution of this system of nonlinear coupled equations is straightforward. Fig. \[condensation\_fintite\_q\]a shows the population $n_{q} = |\psi_{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 = \pi$ is populated on the time scale of $\Gamma^{-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. \[condensation\_fintite\_q\]a) throughout the time-evolution. This contrasts sharply with typical relaxation phenomena in the Brillouin zone (see e.g. Ref. [@HaugJauhoBook]) 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 system-bath coupling ${\cal L}_{\kappa}$.

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 $\psi_{0}$ and $\psi_{\pi}$ of the modes $q = 0$ and $q =
\pi$, respectively, are assumed finite. The EOM (\[eq:boltzmanncoherent\]) then reduce to $$\begin{array}{l}
\partial_{t} \psi_{0} =  - (\Gamma - 4 \kappa ~ \frac{1}{L} |\psi_{\pi}|^{2}) \psi_{0}, \\
\partial_{t} \psi_{\pi} = - i \Omega_{\pi} - [ \Gamma + 4 \kappa ~ (\frac{1}{L}|\psi_{0}|^{2} + 1/2) ] \psi_{\pi}.
\end{array}$$ The structure of these equations is exactly the same as Eq. (\[meanfieldEOM2cavities\]), with the $|a|^2/2$ and $|s|^2/2$ replaced by the photon densities $|\psi_\pi|^2/L$ and $|\psi_0|^2/L$, respectively. The analytical discussion of these coupled equations follows the same path as that from Eq. (\[meanfieldEOM2cavities\]) to Eq. (\[semiclassicalresults\]). The threshold, above which photons start to condense into the $q=0$ mode, is $\Omega_{\pi,{\rm crit}}= \sqrt{L/2}\times \sqrt{\Gamma (\Gamma + 2 \kappa)^2 /
  2 \kappa}$. Note that compared to the results presented in the previous section, a factor $2\kappa$, instead of $\kappa$, 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. \[condensation\_fintite\_q\]b 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, $\langle
\hat{c}_{q}^{\dag} \hat{c}_{q'} \rangle$. 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 $|\langle \hat{c}_{q}^{\dag} \hat{c}_{q'}
\rangle| \simeq \sqrt{n_{q}} \sqrt{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 (\[effectivecavity\]). 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}=\langle
\hat{c}_{q}^{\dag} \hat{c}_{q} \rangle$ are $$\begin{aligned}
\label{cqdagcq}
\partial_{t} \langle \hat{c}_{q}^{\dag} \hat{c}_{q} \rangle &=& -2 [ \Gamma + \kappa (1 - \cos{q})] \langle \hat{c}^{\dag}_{q} \hat{c}_{q} \rangle + i \Omega_{q}^{\ast} \langle \hat{c}_{q} \rangle - i \Omega_{q} \langle \hat{c}_{q}^{\dag} \rangle \nonumber\\ &+& \frac{8 \kappa}{L}  \sum_{q_1} \cos^2\frac{q}{2}  \sin^2\frac{q_1}{2} \langle \hat{c}^{\dag}_{q_1} \hat{c}_{q_1} \rangle + {\cal S}(q),\end{aligned}$$ where the term ${\cal S}(q)$ contains fourth-order correlation functions, which are not diagonal in momentum space, and reads $$\mathcal{S}(q)=   \frac{\kappa}{L} \sum_{q_{1}q_{2}}Ê \left[\cos(q)-\cos(q_1+q_2-q)\right]  \langle c_{q_2+q_1-q}^\dag  c_{q}^\dag c_{q_1}  c_{q_2} \rangle +  {\rm H.c.}$$ 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 $$\label{HFfactorization}
\langle \hat{c}^{\dag}_{q_{1}} \hat{c}^{\dag}_{q_{2}} \hat{c}_{q_{3}} \hat{c}_{q_{4}} \rangle
\simeq (\delta_{q_{1},q_{4}}\delta_{q_{2},q_{3}} + \delta_{q_{1},q_{3}}\delta_{q_{2},q_{4}}) \langle \hat{c}^{\dag}_{q_{1}} \hat{c}_{q_{1}} \rangle \langle \hat{c}^{\dag}_{q_{2}} \hat{c}_{q_{2}} \rangle.$$ This prescription is in fact a Hartree-Fock approximation [@FetterWaleckaBook], 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 $$\mathcal{S}(q)\simeq    \frac{4 \kappa}{L} \sum_{q_{1}}Ê \left[\cos(q)-\cos(q_1)\right]  \langle c_{q_1}^\dag  c_{q_1}\rangle \langle c_q^\dag c_{q}  \rangle.$$ A self-consistent solution is to assume that $\langle \hat{c}_{q} \rangle = 0$ for $q \neq \pi$. The EOM (\[cqdagcq\]) become then a set of coupled, nonlinear, Boltzmann-like equations [@HaugJauhoBook] for the populations only. However, on the edge of the Brillouin zone, the coherent pumping $\Omega_{\pi}$ 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 (\[eq:boltzmanncoherent\]), with a comparable numerical effort. Typical profiles of the steady-state distributions are shown in Fig. \[figurenq\]. 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.

![\[figurenq\] Results of the numerical integration of the Boltzmann-like equations for the populations $n_{q}$ on a discretized Brillouin zone with $L = 201$ points, until $t_{\rm max} \, \Gamma = 10.0$, for $\kappa = 5.0 \, \Gamma$, $\Omega_{\pi} = 3.5 \sqrt{L}\Gamma$ (solid line in (a) and circles in (b)) and $\kappa = 1.0$, $\Omega_{\pi} = 2.0 \sqrt{L}\Gamma$ (dotted line in (a) and crosses in (b)). Two Bose-Einstein equilibrium distributions are plotted in (b) (solid lines) for comparison. ](fig07a "fig:"){width="0.48\linewidth"} ![\[figurenq\] Results of the numerical integration of the Boltzmann-like equations for the populations $n_{q}$ on a discretized Brillouin zone with $L = 201$ points, until $t_{\rm max} \, \Gamma = 10.0$, for $\kappa = 5.0 \, \Gamma$, $\Omega_{\pi} = 3.5 \sqrt{L}\Gamma$ (solid line in (a) and circles in (b)) and $\kappa = 1.0$, $\Omega_{\pi} = 2.0 \sqrt{L}\Gamma$ (dotted line in (a) and crosses in (b)). Two Bose-Einstein equilibrium distributions are plotted in (b) (solid lines) for comparison. ](fig07b "fig:"){width="0.48\linewidth"}

To discuss the width of the central peak, which is related to the correlation length $\xi$ of the condensate (which we give here in units of the lattice spacing), we can restrict our analysis to the dynamics of the $\pi$ mode and to the populations of modes with $q \approx 0$. For the coherence and population of the $q=\pi$ mode we obtain $$\label{picoherence}
\partial_t \langle \hat{c}_\pi \rangle \simeq  - (i\delta_{\rm d}+ \Gamma +2 \kappa+4\kappa {\cal N}_0)  \langle \hat{c}_\pi \rangle -i \Omega_\pi,$$ and $$\partial_t \langle \hat{c}^\dag_\pi \hat{c}_\pi \rangle \simeq  -2\left[\Gamma+2\kappa +4 \kappa {\cal N}_0 \right] \langle \hat{c}^{\dag}_{\pi} \hat{c}_{\pi} \rangle + i \Omega_{\pi}^{\ast} \langle \hat{c}_{\pi} \rangle - i \Omega_{\pi} \langle \hat{c}_{\pi}^{\dag}\rangle.$$ Here we have defined $${\cal N}_0 \equiv \frac{1}{L}\sum'_{q} n_q,$$ with $\sum'_q$ excluding $q=\pi$, and used the Hartree-Fock decoupling approximation, which for three operators reads $\langle
\hat{c}^{\dag}_{q_{1}+q_{2}-q} \hat{c}_{q_{1}} \hat{c}_{q_{2}} \rangle \simeq
n_{q_{1}} \delta(q_{2}-q) \langle \hat{c}_{q} \rangle + n_{q_{2}}
\delta(q_{1}-q) \langle \hat{c}_{q} \rangle$. For $\delta_{\rm d}=0$, the normalized steady-state population ${\cal N}_{\pi} = \langle \hat{c}^\dag_\pi \hat{c}_\pi \rangle/ L$ of the driven mode is then given by $$\label{pidensity}
{\cal N}_{\pi} =   \frac{|\Omega_\pi|^2/L} { \left[ \Gamma +2 \kappa(1+2{\cal N}_0)\right]^2 } .$$ For the low-momentum modes we find $$\partial_t n_q  \simeq   -\left[2\Gamma-  8\kappa{\cal N}_{\pi}+   \kappa(1+2{\cal N}_0+2{\cal N}_{\pi})q^2 \right] n_q+ 8 \kappa(1-q^2/4) {\cal N}_{\pi},$$ and the steady-state solution can be approximately written in the form $$\label{HFsol}
n_{q} \simeq  \frac{n_0}{1 + \xi^{2} q^{2}},$$ where $$\begin{aligned}
\label{n0xi}
n_0= \frac{8 \kappa {\cal N}_{\pi}}{2 \Gamma - 8 \kappa {\cal N}_{\pi}},\qquad \xi^{2} = \frac{\kappa(1+2{\cal N}_0+2{\cal N}_{\pi})}{2 \Gamma - 8 \kappa {\cal N}_{\pi}}+\frac{1}{4}.\end{aligned}$$

Equation (\[HFsol\]) is a Lorentzian profile with broadening $1/\xi$. This scale plays the role of a correlation length since, in the translationally-invariant case, $$\langle \hat{c}^{\dag}_{\ell} \hat{c}_{\ell'} \rangle = \sum_{q,q'} e^{-i q \ell} e^{-i q' \ell'} \langle \hat{c}^{\dag}_{q} \hat{c}_{q'} \rangle \simeq \sum_{q} e^{-i q (\ell - \ell')} \langle \hat{c}^{\dag}_{q} \hat{c}_{q} \rangle \propto e^{-|\ell -\ell'| / \xi}.$$ 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 low-momentum modes are populated by photons which are scattered incoherently from the $q=\pi$ mode.

Discussion
----------

For the Lorentzian low-momentum distribution (\[HFsol\]) we obtain ${\cal N}_0\simeq
n_0/\xi$, and we can interpret the quantity $L\times {\cal N}_0$ as the total number of condensed photons, while the population $L \times {\cal N}_{\pi}$ serves as finite photon reservoir. In steady state, the values of ${\cal N}_{\pi}$, $n_0$ and $\xi$ can be determined from the nonlinear Eqs. (\[pidensity\]) and (\[n0xi\]). The stability of the low-momentum modes requires that ${\cal N}_{\pi} <
\Gamma/4\kappa$, but for a driving strength above $\Omega_{\pi,{\rm crit}}$ the value of ${\cal N}_{\pi}$ will approach this limit. Therefore, in this regime $$\label{Nxi}
{\cal N}_0=\frac{n_0}{\xi}\simeq  \frac{(2\kappa+\Gamma)}{4\kappa} \left(\tilde \Omega_\pi-1\right), \qquad \xi \simeq \frac{\kappa {\cal N}_0 }{2\Gamma} \left( 1+ \frac{\Gamma}{2\kappa} +2 {\cal N}_0 \right),$$ where $\tilde \Omega_\pi=|\Omega_\pi|/\Omega_{\pi,{\rm crit}}>1$. We see that above the critical driving strength the number of photons in the low-momentum modes increases linearly with $\tilde \Omega_\pi$. At the same time the correlation length increases $\sim \tilde \Omega_\pi^2$, indicating an even stronger tendency for photons to occupy the $q=0$ mode. Note that while in the classical high-photon limit, $\Gamma\gg\kappa$, we obtain the ratio ${\cal N}_0/{\cal N}_{\pi} \simeq
(\tilde \Omega_\pi-1)$, in the opposite case, $\kappa\gg \Gamma$, we get an additional enhancement factor, ${\cal N}_0/{\cal N}_{\pi} \simeq (2\kappa/\Gamma)(\tilde
\Omega_\pi-1)$. Therefore, in this regime, very pure, dilute photon condensates with ${\cal N}_0\approx 1$ and $\xi\simeq (\kappa/4\Gamma) (\tilde
\Omega_\pi-1)\tilde \Omega_\pi\gg 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 non-equilibrium open system studied here to a thermal equilibrium gas of massive bosons with finite chemical potential. To do so, we compare the distribution (\[HFsol\]) with a Bose-Einstein distribution $n_q^{\rm
  BE}=1/(\exp\{(\varepsilon_q-\mu)/k_B T\} -1)$ (c.f. Fig. \[figurenq\]b), where $\varepsilon_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 $$\begin{aligned}
 \label{Tmu}
\frac{k_B T}{J} = \frac{1+n_0}{\xi^2}, \qquad
\frac{\mu}{J} = -2 - \frac{1+n_0}{\xi^2} \ln \left( 1+ \frac{1}{n_0} \right).\end{aligned}$$

Note that the temperature $T$ defined in this way is different from $T_{\rm
  eff}$ associated to a finite $\kappa'$, and arises purely from the competition between photon losses and equilibration processes in a driven system. Above threshold, where $n_0,\xi \gg1$, we obtain the scaling $k_BT/J\simeq 4/\tilde
\Omega$ and $k_BT/J\simeq 2\Gamma/(\kappa\tilde \Omega) $ for the limits $\Gamma\gg\kappa$ and $\Gamma\ll \kappa$, respectively. This confirms the conclusion from above that, for $\kappa/\Gamma\rightarrow \infty$, arbitrarily pure photonic condensates can be prepared with our scheme.

As a final remark, we point out that the detuning $\delta_{\rm d}$ of the driving field, which enters as an energy offset in the cavity Hamiltonian (\[Hamiltonianq\]), affects the effective chemical potential $\mu$ only indirectly via a modification of $\Omega_{\pi,{\rm crit}}$. In particular, and in strong contrast to equilibrium systems, the sign of $\delta_{\rm 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. (\[Tmu\]) and (\[Nxi\]). 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 {#ConclusionsSec}
=======================

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 non-trivial 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 non-standard dissipation processes, as well as to study open, non-equilibrium quantum many-body systems.

Acknowledgments {#acknowledgments .unnumbered}
===============

The authors would like to thank P. Zoller for proposing the initial direction of this project and further insightful discussions. We also thank F. Marquardt for valuable feedback on this work. This work was supported by the EU network AQUTE and the Austrian Science Fund (FWF) through SFB FOQUS and the START grants Y 591-N16 (P.R.) and Y 581-N16 (S.D.).

Supplement to Section \[Results\] {#AppendixB}
=================================

In this Appendix we extend Section \[Results\]. Specifically, we study the case of incoherent pumping of the cavities and analyze the stability of the steady-state solutions found in that section.

*Incoherent pumping.* It is important to notice that similar condensation effects to those shown in Section \[Results\] can be achieved under incoherent pumping of the cavity modes. Neglecting the term $\Omega \left(
  \hat{c}_a^{\dagger}+\hat{c}_a \right)$ in the model, and including $$\begin{aligned}
{\cal L}_{\rm pump} = p {\cal D}[\hat{c}_s^{\dagger}] + p {\cal D}[\hat{c}_a^{\dagger}],\end{aligned}$$ we arrive at the EOM $$\begin{array}{l} \label{nEOM2cavities}
\partial_{t} \langle \hat{n}_{\rm s} \rangle = 2 p (\langle \hat{n}_{\rm s} \rangle + 1) - 2 \Gamma \langle \hat{n}_{\rm s} \rangle + 2 \kappa \langle \hat{n}_{\rm a} ( \hat{n}_{\rm s}  + 1) \rangle. \\
\partial_{t} \langle \hat{n}_{\rm a} \rangle = 2 p (\langle \hat{n}_{\rm a} \rangle + 1) - 2 \Gamma \langle \hat{n}_{\rm a} \rangle - 2 \kappa \langle \hat{n}_{\rm a} ( \hat{n}_{\rm s}  + 1) \rangle.
\end{array}$$ Solving for the steady-state solution of these equations within the semiclassical approximation $\langle\hat{n}_s\hat{n}_a\rangle \approx
\langle\hat{n}_s\rangle \langle \hat{n}_a \rangle$, we find the condition $\frac{p}{\Gamma} = \frac{n}{n+2}$, being $n\equiv \langle \hat{n}_s
\rangle+\langle \hat{n}_a\rangle$ the total number of particles in the steady state. Substituting then $\langle
\hat{n}_a\rangle=n-\langle\hat{n}_s\rangle=\frac{2p}{\Gamma-p}-\langle\hat{n}_s\rangle$ in the first steady-state equation, for $p\neq\Gamma$, we find a smooth crossover between an equal population of symmetric and antisymmetric modes for $\kappa\ll p,\Gamma$, and condensation to the symmetric mode for $\kappa\gg
p,\Gamma$ (see Fig. \[incoherent\_pumping\]a). Interestingly, for $p\to\Gamma$, however, all photons occupy the symmetric mode for all values of $\kappa\neq0$. In the latter case, heating (incoherent pumping) compensates cooling, and only the engineered dissipative mechanism survives.

![\[incoherent\_pumping\] a) Fraction of photons in the symmetric mode, $n_s/n\equiv\langle\hat{n}_s\rangle/n$ as a function of the dissipative rate $\kappa$, for incoherent pumping $p/\Gamma=0.5$. b) Threshold given by Eq. (\[threshold\]), below which (\[solbelowthreshold\]) is the only stable solution and above which (\[solabovethreshold\]) is the only stable solution.](fig08a "fig:"){width="0.48\linewidth"} ![\[incoherent\_pumping\] a) Fraction of photons in the symmetric mode, $n_s/n\equiv\langle\hat{n}_s\rangle/n$ as a function of the dissipative rate $\kappa$, for incoherent pumping $p/\Gamma=0.5$. b) Threshold given by Eq. (\[threshold\]), below which (\[solbelowthreshold\]) is the only stable solution and above which (\[solabovethreshold\]) is the only stable solution.](fig08b "fig:"){width="0.48\linewidth"}

*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 (\[solbelowthreshold\]) and (\[solabovethreshold\]). To this end we linearize the EOM (\[meanfieldEOM2cavities\]) as $s\to s_0 + \delta s$, $a\to
a_0 + \delta a$ (similarly for the conjugate fields), being $s_0$, $a_0$ the steady-state solutions and $\delta s$, $\delta a$ representing fluctuations. This gives $\frac{d}{dt} (\delta s, \delta a, \delta s^* \delta
a^*)_T = J (\delta s, \delta a, \delta s^* \delta a^*)_T$ with the $4\times
4$ Jacobian matrix $J$: $$\left ( \begin{array}{cccc}
-(\Gamma-\kappa|a_0|^2) & \kappa |s_0| a_0^* & 0 & \kappa |s_0| a_0 \\
-\kappa |s_0| a_0 & i\delta_c- [\Gamma+\kappa(|s_0|^2+1)] & -\kappa |s_0| a_0 & 0 \\
0 & \kappa |s_0| a_0^* & -(\Gamma-\kappa |a_0|^2) & \kappa |s_0| a_0 \\
-\kappa |s_0| a_0^* & 0 & - \kappa |s_0| a_0^* & -i\delta_c- [\Gamma+\kappa(|s_0|^2)]
\end{array} \right ),$$ 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 steady-state solutions presented in Section \[Results\] are stable in the whole range of parameters. In Fig. \[incoherent\_pumping\]b we show the threshold $\Omega_{\rm crit}$ given by Eq. (\[threshold\]), below which (\[solbelowthreshold\]) is the only stable solution and above which (\[solabovethreshold\]) is the only stable solution.

References {#references .unnumbered}
==========

[99]{}

Deng H, Haug H and Yamamoto Y 2010 [*Rev. Mod. Phys.*]{} [**82**]{} 1489

Deng H [*et al*]{} 2002 [*Science*]{} [**298**]{} 199

Kasprzak [*et al*]{} 2006 [*Nature*]{} [**443**]{} 409

Balili R [*et al*]{} 2007 [*Nature*]{} [**316**]{} 1007

Amo A [*et al*]{} 2009 [*Nature*]{} [**457**]{} 291

Bloch I, Dalibard J and Zwerger W 2008 [*Rev. Mod. Phys.*]{} [**80**]{} 885

Klaers J, Vewinger F and Weitz M 2010 [*Nature Phys.*]{} [**6**]{} 512

Klaers [*et al*]{} 2010 Nature [**468**]{} 545

Klaers J [*et al*]{} 2012 [*arXiv:1201.0444*]{}

Snoke D W 2003 [*Phys. Stat. Sol (b)*]{} [**238**]{} 389

Szymanska M H and Littlewood P B 2002 [*Solid State Comm.*]{} [**124**]{} 103

Szymanska M H, Littlewood P B and Simons B D 2003 [*Phys. Rev. A*]{} [**68**]{} 013818

Hartmann M J, Brandao F G S L and Plenio M P 2008 [*Laser Photonics Rev.*]{} [**2**]{} 527

Tomadin A and Fazio R 2010 [*J. Opt. Soc. Am. B*]{} [**27**]{} A130

Hartmann M J, Brandao F G S L and Plenio M B 2006 [*Nature Phys.*]{} [**2**]{} 849

Greentree A D, Tahan C, Cole J H and Hollenberg C L 2006 [*Nature Phys.*]{} [**2**]{} 856

Angelakis D G, Santos M F and Bose S 2007 [*Phys. Rev. A*]{} [**76**]{} 031805(R)

Rossini D and Fazio R 2007 [*Phys. Rev. Lett.*]{} [**99**]{} 186401

Na N [*et al*]{} 2008 [*Phys. Rev. A*]{} [**77**]{} 031803(R)

Fleischhauer M, Otterbach J and Unanyan R G 2008 [*Phys. Rev. Lett.*]{} [**101**]{} 163601

Koch J and Le Hur K 2009 [*Phys. Rev. A*]{} [**80**]{} 023811

Schmidt S [*et al*]{} 2010 [*Phys. Rev. B*]{} [**82**]{} 100507

Raimond J M, Brune M and Haroche S 2001 [*Rev. Mod. Phys.*]{} [**73**]{} 565

Imamoglu A 1997 [*Phys. Rev. Lett.*]{} [**79**]{} 1467

Birnbaum K M [*et al*]{} 2005 [*Nature*]{} [**436**]{} 87

Hennessy K [*et al*]{} 2007 [*Nature*]{} [**445**]{} 896

Faraon A [*et al*]{} 2008 [*Nature Phys.*]{} [**4**]{} 859

Rebić S, Twamley J and Milburn G J 2009 [*Phys. Rev. Lett.*]{} [**103**]{}, 150503

Cho J, Angelakis D G and Bose S 2008 [*Phys. Rev. Lett.*]{} [**101**]{} 246809

Koch J [*et al*]{} 2010 [*Phys. Rev. A*]{} [**82**]{} 043811

Hafezi M [*et al*]{} 2011 [*Nature Phys.*]{} [**7**]{} 907

Umucalilar R O and Carusotto I 2011 [*arXiv:1110.6524*]{}

Chang D E [*et al*]{} 2008 [*Nature Phys.*]{} [**4**]{} 884

Carusotto I [*et al*]{} 2009 [*Phys. Rev. Lett.*]{} [**103**]{} 033601

Lepert G [*et al*]{} 2011 [*New Journ. Phys.*]{} [**13**]{} 113002

Wallraff A [*et al*]{} 2004 [*Nature*]{} [**431**]{} 162

Chiorescu I [*et al*]{} 2004 [*Nature*]{} [**431**]{} 159

Blais A [*et al*]{} 2004 [*Phys. Rev. A*]{} [**69**]{} 062320

Schoelkopf R J and Girvin S M 2008 [*Nature*]{} [**451**]{} 664

Clarke J and Wilhelm F K 2008 [*Nature*]{} [**453**]{} 1031

Mariantoni M [*et al*]{} 2011 [*Nature Phys.*]{} [**7**]{} 287

Tomadin A [*et al*]{} 2010 [*Phys. Rev. A*]{} [**81**]{} 061801(R)

Pippan P, Evertz H G and Hohenadler M 2009 [*Phys. Rev. A*]{} [**80**]{} 033612

Knap M, Arrigoni E and von der Linden W 2010 [*Phys. Rev. B*]{} [**81**]{} 104303

Leib M and Hartmann M 2010 [*New Jour. Phys.*]{} [**12**]{} 093031

Knap M [*et al*]{} 2011 [*Phys. Rev. A*]{} [**83**]{} 023821

Poyatos J F, Cirac J I and Zoller 1996 P [*Phys. Rev. Lett.*]{} [**77**]{} 4728

Diehl S [*et al*]{} 2008 [*Nature Phys.*]{} [**4**]{} 878

Verstraete F, Wolf M M and Cirac J I 2009 [*Nature Phys.*]{} [**5**]{} 633

Weimer H 2010 [*et al*]{} [*Nature Phys.*]{} [**6**]{} 382

Kastoryano M J, Reiter F and Sørensen A S 2011 [*Phys. Rev. Lett.*]{} [**106**]{} 090502

Muschik C A, Polzik E S and Cirac J I 2011 [*Phys. Rev A*]{} [**83**]{} 052312

Barreiro J T [*et al*]{} 2011 [*Nature*]{} [**470**]{} 486

Krauter H [*et al*]{} 2011 [*Phys. Rev. Lett.*]{} [**107**]{} 080503

Diehl S, Tomadin A, Micheli A, Fazio R and Zoller P 2010 [*Phys. Rev. Lett.*]{} [**105**]{} 015702

Tomadin A, Diehl S and Zoller P 2011 [*Phys. Rev. A*]{} [**83**]{} 013611

Devoret M 1997 [*in: QuantumFluctuations, edited by S.Reynaud, E. Giacobino and J. Zinn-Justin*]{} (Amsterdam: Elsevier)

Makhlin Y, Schön G and Shnirman A 2001 [*Rev. Mod. Phys.*]{} [**73**]{} 357

Breuer H P and Petruccione F 2002 [*The Theory of Open Quantum Systems*]{} (New York: Oxford)

Gardiner C W and Zoller P 2000 [*Quantum Noise*]{} (Berlin: Springer)

Watanabe G and Mäkelä H 2011 [*arXiv:1101.4845*]{}

Haken H 1985 [*Light – Laser Light Dynamics*]{} (Amsterdam: North-Holland)

Astafiev O [*et al*]{} 2007 [*Nature*]{} [**449**]{} 588

Rodrigues D A, Imbers J and Armour A D 2007 [*Phys. Rev. Lett.*]{} [**98**]{} 067204

You J Q [*et al*]{} 2007 [*Phys. Rev. B*]{} [**75**]{} 104516

Hauss J [*et al*]{} 2008 [*Phys. Rev. Lett.*]{} [**100**]{} 037003

Scully M O 1999 [*Phys. Rev. Lett.*]{} [**82**]{} 3927

Kocharovsky V V [*et al*]{} 2000 [*Phys. Rev. Lett.*]{} [**61**]{} 023609

Haug H and Jauho A-P 1997 [*Quantum Kinetics in Transport and Optics of Semiconductors*]{} (New York: Springer)

Fetter A L and Walecka J D 2003 [*Quantum Theory of Many-Particle Systems*]{} (New York: Dover)

[^1]: The explicit derivation using circuit theory is omitted here. For a comprehensive review on the topic c.f. [@BHWGS2004; @D97; @MSS2001].