---
abstract: 'We propose and study the use of photon-mediated interactions for the generation of long-range steady-state entanglement between $N$ atoms. Through the judicious use of coherent drives and the placement of the atoms in a network of Cavity QED systems, a balance between their unitary and dissipative dynamics can be precisely engineered to stabilize a long-range correlated state of qubits in the steady state. We discuss the general theory behind such a scheme, and present an example of how it can be used to drive a register of $N$ atoms to a generalized W state and how the entanglement can be sustained indefinitely. The achievable steady-state fidelities for entanglement and its scaling with the number of qubits are discussed for presently existing superconducting quantum circuits. While the protocol is primarily discussed for a superconducting circuit architecture, it is ideally realized in any Cavity QED platform that permits controllable delivery of coherent electromagnetic radiation to specified locations.'
author:
- Camille Aron
- Manas Kulkarni
- 'Hakan E. Türeci'
bibliography:
- 'SpinChain.bib'
title: |
    Photon-mediated interactions: a scalable tool\
    to create and sustain entangled states of $N$ atoms
---

Introduction {#sec:intro}
============

Photon-mediated interactions are ubiquitous in nature. While the traditional formulation of quantum electrodynamics places equal emphasis on fields and sources, it is possible to take a point of view where the electromagnetic degrees of freedom are integrated out to reach an effective nonlocal field theory for matter only [@feynman_mathematical_1950]. A particularly beautiful example of this point of view, most closely related to the phenomena investigated here, is Schwinger’s formulation of the Casimir effect [@schwinger_casimir_1975]. Here, the electromagnetic degrees of freedom are integrated out and result in a photon-mediated retarded and nonlocal interaction between two conducting surfaces.

A classic example where photon-mediated interactions are at play is the superradiance (and subradiance) of a cluster of dipoles, first pointed out by Dicke [@dicke_coherence_1954]. Here, photon-mediated interactions ultimately lead to generation of transient coherence between dipoles [@gross_superradiance_1982; @brandes_coherent_2005], which results in the emission of a powerful pulse whose intensity scales with $N^2$, where $N$ is the number of dipoles within a volume $V \sim \lambda^3$ ($\lambda$ is the wavelength of radiation). In free space, however, such interactions decay fast with interdipole distance, and it is challenging to engineer such interactions in a controlled way [@devoe_observation_1996; @eschner_light_2001]. Photon-mediated interactions are also at play in the self-organization transition of optically driven cold atoms in a cavity [@ritsch_cold_2013]. In these systems, cavity-mediated long-range interactions between atoms, tunable by the drive strength, lead to softening of a motional excitation mode recently observed in experiments [@baumann_dicke_2010; @mottl_roton-type_2012]. Certain aspects of the underlying critical behavior of this non-equilibrium phase transition can be described through photon-mediated interactions between the atoms constituting the condensate [@brennecke_real-time_2013; @strack_dicke_2011; @kulkarni_cavity-mediated_2013; @konya_damping_2014].

In recent years, we have seen the first attempts to use such photon-mediated interactions to generate strong coupling and possibly entanglement between artificial atoms in solid-state cavity QED systems. These approaches capitalize on strong light-matter interactions that can be generated in confined geometries such as resonators [@majer_coupling_2007; @filipp_multimode_2011] and waveguides [@loo_photon-mediated_2013]. In particular, Ref. [@loo_photon-mediated_2013] demonstrated coherent exchange interactions between two superconducting qubits separated by as much as a full wavelength (in that case $\lambda \sim 18.6$ mm) in an open quasi-1D transmission line. More recently, superradiance of two artificial atoms was observed and characterized in a controlled setting in a superconducting quantum circuit [@mlynek_observation_2014].

![(color online) (a) One-dimensional array of cavity-qubit systems coupled by photon exchange and subject to one or several ac microwave drives, cavity decay at rate $\kappa$, qubit relaxation $\gamma$ and pure dephasing $\gamma_\phi$. (b) Implementation with superconducting transmon qubits embedded in interconnected microwave cavities and driven by external continuous-wave generators.[]{data-label="figschem1d"}](many1d.eps "fig:"){width="8.5cm"} ![(color online) (a) One-dimensional array of cavity-qubit systems coupled by photon exchange and subject to one or several ac microwave drives, cavity decay at rate $\kappa$, qubit relaxation $\gamma$ and pure dephasing $\gamma_\phi$. (b) Implementation with superconducting transmon qubits embedded in interconnected microwave cavities and driven by external continuous-wave generators.[]{data-label="figschem1d"}](irf.eps "fig:"){width="7.0cm"}

The goal of this paper is to show how photon-mediated interactions between qubits embedded in an engineered electromagnetic environment can be harnessed to controllably generate a certain large-scale entangled state of $N$ qubits [*in the steady state*]{}. In the present work, the role of photons is twofold. First, they mediate a coherent coupling between distant qubits. Second, they provide a controllable dissipative mechanism which can be used to stabilize a long-range correlated many-body state of the qubits. By driving the system at a suitable frequency, one can achieve a transition that produces the desired many-body state while dumping energy into one of the electromagnetic modes. The dissipative mechanisms are key to make the scheme steady-state, in contrast to better-known alternatives, such as those based on Rabi cycling between a ground state and a desired excited state. We show how the balance of the unitary and dissipative contributions can be precisely tuned by the placement of qubits in an engineered photonic environment, and a set of coherent electromagnetic drives with specified amplitudes and frequencies.

Engineered dissipative dynamics has recently been employed in superconducting circuits to cool a single qubit to a desired state on the Bloch sphere [@murch_cavity-assisted_2012], and two qubits that reside in a single cavity to a Bell state [@shankar_autonomously_2013]. Furthermore, a number of recent theoretical works have focused on generation of high-fidelity steady-state entanglement between two superconducting qubits [@leghtas_stabilizing_2013; @reiter_steady-state_2013; @huang_generation_2013; @aron_steady-state_2014]. We present the general theory underlying these phenomena, deriving and solving the dynamics of $N$ qubits residing in an arbitrary open electromagnetic environment, subject to coherent driving and losses. [The underlying new principle is based on using collective electromagnetic (EM) modes of a large structure (such as Bloch modes of a lattice of cavities) to dissipatively stabilize a collective state of spins. By “scalability” we refer to the fact that, for the specific protocol worked out here for W states, that the same protocol used for $N$ qubits can be used for $N+1$ qubits. Generally it is to be expected that the fidelity of stabilization degrades as $N$ is increased because of spectral crowding, which we discuss for the particular protocol we investigate here, and provide an estimate for the maximum number of qubits that can be stabilized reliably. Interestingly, spatial symmetries in the system (such as discrete translational invariance) can be used along with a [*spatially*]{} modulated coherent drive, to greatly enhance the spectral resolution. In one of the protocols we propose, this allows us to resolve W states in a way that the state selectivity is limited by only the finesse of the EM system, instead of the mean level spacing of the collective spin states.]{}

In Sec. \[sec:1dmodel\], we lay out a particular architecture that we have in mind for a proof-of-principle demonstration: we consider of model of $N$ qubits residing in a one-dimensional array of cavities. We present the different layers of approximations that allow us to have an analytic handle over the problem. First, we show how cavity photons mediate effective qubit-qubit interactions through the collective electromagnetic modes of the array. In Sec. \[sec:cooling\], we show how the [*full Liouvillian*]{} describing the evolution of the reduced density matrix of the $N$ qubits can be engineered to drive the qubit subsystem to generalized W states, and the entanglement sustained as long as the drives are on. The fabrication tolerances of such a cavity array has recently been studied experimentally [@underwood_low-disorder_2012]. Thus, the controlled fabrication of such arrays is well within reach of the current superconducting circuit technology both in coplanar and 3D configurations [@underwood_thesis_2015]. We present the fidelities that can be expected in recently fabricated systems and analyze the fault tolerance of the method to phase and amplitude noise of the drive parameters, as well as the nonuniformity of qubit and cavity parameters. Finally, in Sec. \[sec:ph-mediated\], we generalize our scheme to arbitrary arrangements of qubits coupled to an engineered photonic backbone, going beyond the tight-binding approximation for the EM system used in previous sections.

One-dimensional lattice model {#sec:1dmodel}
=============================

Consider a one-dimensional array of $N$ identical microwave cavities with nominally equal frequencies $\omega_{{\rm c}}$, coupled to each other capacitively described by a nearest-neighbor tunneling matrix element $J$. Each cavity houses a superconducting qubit with splitting $\omega_{{\rm q}}$. We shall consider both the cases of open (*i.e.* non-periodic) and periodic boundary conditions (identifying site $N$ with site 0). Later, we will relax our assumptions on identical cavity resonance frequencies and consider the most general case, showing that the general approach to the dissipative stabilization of a generalized W state stands.

Furthermore, each cavity shall be driven by a coherent monochromatic microwave source with frequency $\omega_{{\rm d}}$ and a site-dependent amplitude $\epsilon^{{{\rm d}}}_i$ and phase $\Phi_i$; see Fig. \[figschem1d\](b) for a 3D superconducting circuit architecture of the system Fig. \[figschem1d\](a). We work in a regime where $\omega_{{\rm c}}$, $\omega_{{\rm q}}$ and $\omega_{{\rm d}}$ are mutually far detuned from each other (typically on the order of GHz). The starting Hamiltonian is that of a one-dimensional driven Jaynes-Cummings lattice model studied before in Refs. [@knap_emission_2011; @nissen_nonequilibrium_2012; @grujic_non-equilibrium_2012; @hur_many-body_2015], $$\begin{aligned}
 \label{eq:H}
 H(t) = H_\sigma + H_{\sigma a} + H_a(t) \,,
 \end{aligned}$$ where $H_\sigma$, $H_{\sigma a}$, and $H_a(t)$ are respectively the qubit, the Jaynes-Cummings light-matter coupling, and the driven cavity Hamiltonians, $$\begin{aligned}
H_\sigma =& \! \sum_i   \omega_{{\rm q}}\frac{\sigma_i^z}{2},\,    H_{\sigma a} =  g \! \sum_i\! \left[ a_i^\dagger \sigma^-_i + \mbox{H.c.} \right]  \,,
\\
H_a(t) =& \! \sum_i\!  \left[ \omega_{{\rm c}}a_i^\dagger a_i 
   - J ( a_{i}^\dagger a_{i+1} + \mbox{H.c.} ) \right.
      \nonumber \\
  &   \qquad \left. +  2 \epsilon^{{{\rm d}}}_i \cos(\omega_{{\rm d}}t + \Phi_i) \left( a_i + a^\dagger_i \right)\right] \,. \label{eq:Ha}\end{aligned}$$ Here, $i$ runs from site $0$ to $N-1$. The qubits are two-level systems described by the usual Pauli pseudo-spin operators $\sigma_i^{x,y,z}$ and $\sigma_i^\pm \equiv (\sigma^x_i \pm {{\rm i}}\sigma_i^y)/2$. $\epsilon^{{\rm d}}_i$ and $\Phi_i$ are, respectively, the cavity-dependent amplitude and phase of the ac microwave drives. Thermal equilibrium is achieved by setting all drive amplitudes to zero, $\epsilon^{{\rm d}}_i=0$. Without loss of generality, we assume that the detuning between cavity and qubit frequency $\Delta \equiv \omega_{{\rm q}}- \omega_{{\rm c}}> 0$, and $J>0$. We operate in the dispersive regime corresponding to $g/\Delta \sim 10^{-1}$ and at sufficiently weak drive amplitudes to ensure the presence of very few photons in the cavities, so the Schrieffer-Wolff perturbation theory to be applied shortly is well justified. Below, we entirely integrate out the photonic degrees of freedom, resulting in (a) an effective qubit-qubit interaction as discussed on more general grounds in Sec. \[sec:ph-mediated\], (b) local Zeeman fields of the form $\epsilon^{{\rm d}}_i \sigma_i^{x}$ (for $\Phi_i=0$), and (c) a non-unitary evolution characterized by controllable transition rates.

#### Intrinsic dissipation in the system. {#intrinsic-dissipation-in-the-system. .unnumbered}

In addition to the unitary dynamics described by Eq. (\[eq:Ha\]), we assume the individual qubits are coupled to uncontrolled environmental degrees of freedom that give rise to single qubit spin-flip rate ($\gamma$), a single-qubit pure dephasing rate $\gamma_\phi$, and a cavity decay rate $\kappa$.

#### Typical system parameters. {#typical-system-parameters. .unnumbered}

In a recent experiment studying the $N=2$ case of the protocol described here [@ucbpaper] in a 3D superconducting circuit architecture, typical system parameters were: $\omega_{{\rm c}}\simeq 7, \omega_{{\rm q}}\simeq 6, g \simeq J \simeq 10^{-1}, \kappa \simeq 7 \times 10^{-4}, \gamma \simeq 4 \times 10^{-5}$ all in units of $2\pi$ GHz. Below, our analytical approach is performed assuming the hierarchy of energy scales $\Delta \gg g, J \gg \kappa \gg \gamma \gg \gamma_\phi$.

#### Rotating wave approximation. {#rotating-wave-approximation. .unnumbered}

We eliminate the time dependence in $H_a(t)$ by working in the frame rotating at $\omega_{{\rm d}}$ and dropping the non-secular terms. In the rest of the Hamiltonian (\[eq:H\]), this also amounts to replacing $\omega_{{\rm q}}$ by $\Delta_{{\rm q}}\equiv \omega_{{\rm q}}- \omega_{{\rm d}}$ and $\omega_{{\rm c}}$ by $-\Delta_{{\rm c}}\equiv \omega_{{\rm c}}- \omega_{{\rm d}}$. Once expressed in the eigenbasis of the (undriven) coupled cavity system, $H_{a}$ is given by $$\begin{aligned}
 \label{eq:Hak}
  H_{a} = \sum_k (\omega_k -\omega_{{\rm d}}) a^\dagger_k a_k  +  ( \epsilon^{{\rm d}}_k a^\dagger_k + \mbox{H.c.} )\,.\end{aligned}$$ Here, $k$ is the discrete Bloch wavevector, $a_k^\dagger = \sum_j \varphi_{k}^*(j) \, a_j^\dagger$ creates a photon in the $k^{\mathrm{th}}$ mode. The specific discrete values of $k$ and the corresponding mode profiles $\varphi_k(j)$ are to be fixed by the boundary conditions. For periodic boundary conditions, the set of quasi-momenta are $k= 2\pi\, n/N$, with $n = 0,\ldots, N-1$ and $\varphi_{k}(j) = {{\rm e}}^{-{{\rm i}}k j} / \sqrt{N}$. For open boundary conditions, $ k = \pi (n+1)/(N+1)$ and $\varphi_k(j) = \sqrt{2/(N+1)} \sin\left(k(j+1)\right)$. In Eq. (\[eq:Hak\]), we incorporated the phase $\Phi_j$ in $\epsilon^{{\rm d}}_j$, which can henceforth be complex and $\epsilon^{{\rm d}}_k =   \sum_j \varphi_{k}(j) \, \epsilon_j^{{\rm d}}$. The eigenfrequencies are given by the photonic dispersion relation $$\begin{aligned}
\omega_k = \omega_{{\rm c}}- 2 J \cos(k)\,.\end{aligned}$$

![(color online) Engineered qubit many-body spectrum. The non-equilibrium drive and the resulting photon fluctuation bath are used to create a dominant transition from the trivial ground state to a target entangled W state in the 1-excitation manifold. The existence of a non-zero spontaneous decay $\gamma$ is critical to avoid populating higher-excited state manifolds.[]{data-label="fig:mediated"}](mechaN.eps){width="3.5cm"}

#### Effective dissipative $XY$ model. {#effective-dissipative-xymodel. .unnumbered}

We eliminate the light-matter interaction to second-order perturbation theory in $g/\Delta$ by means of a Schrieffer-Wolf (SW) transformation which maps $H \mapsto {{\rm e}}^{X} H  {{\rm e}}^{X^\dagger}$ where $$\begin{aligned}
X \equiv g \sum_k \left[ \frac{a_k \sigma_k^\dagger}{\omega_{{\rm q}}- \omega_k}
-\mathrm{H.c.}
 \right]\,.\end{aligned}$$ Collecting all the terms, we obtain an isotropic $XY$ model subject to a magnetic field for the qubit subsystem ($H_\sigma$), weakly coupled to the photon fluctuations of the EM backbone ($H_{\sigma a}$): $$\begin{aligned}
H_\sigma =& \! \sum_i  
\boldsymbol{h_i} \cdot \frac{\boldsymbol{\sigma}_i}{2} - \frac{J}{2} \left(\frac{g}{\Delta}\right)^2 \left[ \sigma^x_{i} \sigma^x_{i+1} + \sigma^y_{i} \sigma^y_{i+1}\right]\,, \label{eq:Hspinchain} \\
H_{\sigma a} &= \sum_i  \left(\frac{g}{\Delta}\right)^2 \sigma_i^z( \Delta a_i^\dagger a_i + \epsilon^{{\rm d}}_i a_i^\dagger +{\epsilon^{{\rm d}}_i}^* a_i)\,.
\label{eqtransXY}\end{aligned}$$ Here, ${h}^x_i = 2 \mbox{Re}(\epsilon^{{\rm d}}_i) (g/\Delta)$, $h^y_i= -2 \mbox{Im}(\epsilon^{{\rm d}}_i) (g/\Delta)$, $h^z_i  = \Delta_{{\rm q}}+ \Delta (g/\Delta)^2 $. The effective magnetic field $\boldsymbol{h}_i$ is mainly oriented along the $z$ direction but we shall see that, while they break the integrability of the model, the $x$ and $y$ components of the emergent Zeeman field play a crucial role for the scheme below.

We note that the SW transformation is also responsible for subleading corrections to the strength of the dissipative terms (*e.g.* the qubit decay $\gamma$ acquires a contribution from the cavity decay $\kappa$, corresponding to a Purcell contribution, and *vice versa*) the local contributions of which can be simply included by working with renormalized parameters $\kappa$, $\gamma$ and $\gamma_\phi$ at low photon occupations [@boissonneault_dispersive_2009], which is the regime studied here. These effective parameters can be extracted from experimental measurements.

When the drives are off, $\epsilon^{{\rm d}}_i = 0$, the system simply thermalizes with its environment – this is just the physics of blackbody radiation in a coupled cavity system (including $N$ dipoles). Viewed from the perspective of the qubits, while $H_\sigma$ describes a quantum phase transition from a paramagnetic to a ferromagnetic phase when the magnitude of the transverse field is on the order of the nearest-neighbor coupling $J ( g/\Delta)^2$, this regime is never reached for realistic system parameters. For Raman-driven qubits the situation is more interesting, and the phase diagram was recently studied displaying various exotic attractors [@schiro_exotic_2015]. Henceforth, we only work in the experimentally achievable regime $\omega_{{\rm q}}\gg  J(g/\Delta)^2$ and, to leading order, the ground state of $H_\sigma$ is simply the separable state $|\boldsymbol{0} \rangle  \equiv | \! \downarrow \ldots \downarrow \rangle$. The low-energy spectral content of the qubit sector is depicted in Fig. \[fig:mediated\]: the first $N$-qubit excited manifolds are roughly separated by $\Delta_{{\rm q}}$ while the lifting of degeneracy within each manifold is controlled by $J (g/\Delta)^2$

From the point of view of equilibrium statistical mechanics, the many-body system of qubits described by $H_\sigma$ is not interesting because it would thermalize with the radiative reservoir, reaching a steady state $\rho_\sigma (t \rightarrow \infty) \propto {{\rm e}}^{-\beta H_\sigma} \sim | \boldsymbol{0}  \rangle \langle \boldsymbol{0}|$, a collection of uncorrelated qubits in their ground state. We shall see that turning on the drives will change this situation dramatically. A careful choice of drive parameters (even for small amplitudes of drives) will be shown to enable the stabilization of a particular many-body state of qubits in the excited-state manifold.

Stabilization of generalized W states {#sec:cooling}
=====================================

Our goal in this section is to identify non-trivial entangled eigenstates of the spin chain $H_\sigma$ and design a protocol which, starting from the ground state $|\boldsymbol{0}\rangle$ that can be straightforwardly prepared, achieves the stabilization of an interesting excited state of choice. Below we discuss the details of a robust protocol for the stabilization of a generalized W state of qubits with [*minimal*]{} resources.

#### Linearized photon spectrum. {#linearized-photon-spectrum. .unnumbered}

We first address the non-linearities in $H_{\sigma a}$ by decomposing the photonic field into mean-field plus bosonic fluctuations: $$\begin{aligned}
a_k \equiv \bar a_k + d_k \mbox{ with }  \bar a_k = \frac{\epsilon^{{\rm d}}_k}{\omega_{{\rm d}}- \omega_k + {{\rm i}}\kappa/2}\;. \label{eqsmallfluct}\end{aligned}$$ We assume here that all Bloch modes have the same loss $\kappa$. This can be made more precise but it will not qualitatively change the results we present. Neglecting those terms that are quadratic in the fluctuations and that couple to the qubits, the light sector reduces to $$\begin{aligned}
H_a \rightarrow H_d &= \sum_k (\omega_k - \omega_{{\rm d}}) d_k^\dagger d_k \,, \label{eq:bath1} \\
H_{\sigma a} \rightarrow H_{\sigma d} &= \left(\frac{g}{\Delta}\right)^2  \sum_i  \sigma^z_i [( \Delta \bar{a}_i  + \epsilon^{{\rm d}}_i ) d_i^\dagger + \mbox{H.c.} ]\,. \label{eq:bath2}\end{aligned}$$

#### Diagonalization of the matter sector. {#diagonalization-of-the-matter-sector. .unnumbered}

The Hamiltonian Eq. (\[eq:Hspinchain\]) in the presence of non-zero drive terms is non-integrable. We therefore proceed by projecting it into the low-energy sector with a maximum one excitation: $$\begin{aligned}
 \label{eq:Hstrunc}
H_\sigma  =\sum_k E_k |k \rangle \langle k | + \left( \frac{g}{\Delta} \right) \left( \epsilon^{{\rm d}}_k  |k\rangle\langle \boldsymbol{0} | + \mbox{H.c.} \right),\end{aligned}$$ We have set the energy of the ground state $|\boldsymbol{0} \rangle  \equiv | \downarrow \ldots \downarrow \rangle$ to zero ($E_{\boldsymbol{0}} = 0$). Here the states $$|k\rangle = \sum_{i=0}^{N-1}  \varphi_k^* (i) \ |i\rangle\,,$$ with $|i\rangle \equiv |\downarrow_0 \ldots\downarrow_{i-1} \, \uparrow_i \, \downarrow_{i+1} \ldots \downarrow_{N-1}  \rangle$ indicating one excitation located at site $i$, are states that carry a single qubit excitation of quasi-momentum $k$, entangled over the entire chain. These are the eigenstates of the [*undriven*]{} spin chain (*i.e.* for $\epsilon^{{\rm d}}_i =0 \, \forall\, i$) with a dispersion relation $$\begin{aligned}
E_k = \epsilon_k - \omega_{{\rm d}}, \; \epsilon_k = \omega_{{\rm q}}+ \delta\omega_{{\rm q}}- 2 J \left(\frac{g}{\Delta}\right)^2 \! \cos(k)\;.\end{aligned}$$ Here $\delta\omega_{{\rm q}}\simeq  g^2/\Delta$ is the cavity-induced Stark shift. This truncation of the Hamiltonian holds if the higher-excitation manifolds are not significantly occupied during the dynamics. This can be checked a posteriori and we do so. Let us first discuss the case of open boundary conditions for which the absence of translational and space-reversal symmetry generically yields a fully non-degenerate spectrum. The effect of the drive term in Eq. (\[eq:Hstrunc\]), assumed to be small as stated before, can be taken into account through a perturbation theory and yields the following eigenstates of $H_\sigma$ to lowest order in $({g}/{\Delta}) (\epsilon_k^{{\rm d}}/\Delta_{{\rm q}})$: $$\begin{aligned}
|\widetilde{\boldsymbol{0}} \rangle 
 &\simeq
  |\boldsymbol{0} \rangle \!-\! \left(\frac{g}{\Delta}\right) \!  \sum_{k}  \! \frac{\epsilon^{{\rm d}}_k}{\Delta_{{\rm q}}} |k \rangle,\, 
\widetilde{E}_{\boldsymbol{0}} \simeq \!-\!  \left(\frac{g}{\Delta}\right)^2 \! \sum_k \! \frac{|\epsilon^{{\rm d}}_k|^2}{\Delta_{{\rm q}}},\\
|\widetilde{k} \rangle & \simeq
|k \rangle + \left(\frac{g}{\Delta}\right)  \frac{{\epsilon^{{\rm d}}_k}^*}{\Delta_{{\rm q}}} |\boldsymbol{0} \rangle,\,
\widetilde{E}_k \simeq  E_k + \left(\frac{g}{\Delta}\right)^2 \frac{|\epsilon^{{\rm d}}_k|^2}{\Delta_{{\rm q}}}.
\label{eqqbitcollen}\end{aligned}$$ The above corrections to the undriven eigenstates are crucial for the success of the two-photon cooling mechanism presented below.

#### Transition rates. {#transition-rates. .unnumbered}

By virtue of Eq. (\[eqsmallfluct\]), the coupling of the photonic fluctuations (on top of a classical part) to the spin-chain Eq. (\[eqtransXY\]) can be treated in perturbation theory. This permits us to integrate them out arriving at an effective master equation for qubits only. We note that the fluctuations of the collective photon modes of the lattice, described by spectral function per mode $q$ $$\begin{aligned}
\rho_q(\omega) = - \frac{1}{\pi} \mbox{Im } \frac{1}{\omega -\omega_q+{{\rm i}}\kappa/2}\;,\end{aligned}$$ can be easily manipulated by the design of the cavity lattice. Assuming that the photon fluctuations, which couple to the spin degrees of freedom via [Eq. (\[eq:bath2\])]{}, thermalize with the radiative environment which in turn is taken to be at very low temperature, we arrive at the effective master equation for spins only $\rho_\sigma$ at steady state, $\rho_\sigma^{\rm NESS} \equiv \lim\limits_{t\to\infty} \rho_\sigma$: $$\begin{aligned}
\label{eq:master}
\partial_t \rho^{\rm NESS}_\sigma & \! \! =  0 =  -{{\rm i}}\left[ H_\sigma,\rho^{\rm NESS}_\sigma \right]
+ 
\sum_k \Gamma_{\boldsymbol{0}\to k} \mathcal{D}[| \widetilde{k} \rangle \langle \widetilde{\boldsymbol{0}} | ] \rho^{\rm NESS}_\sigma \nonumber \\
& \hspace{-3em} + \gamma \sum_k  \mathcal{D}[| \widetilde{\boldsymbol{0}} \rangle \langle \widetilde{k} | ]  \rho^{\rm NESS}_\sigma
 +  \frac{2\gamma_\phi}{N} \sum_{k\,q} \mathcal{D}[| \widetilde{q} \rangle \langle \widetilde{k} | ] \rho^{\rm NESS}_\sigma.\end{aligned}$$ The derivation of this master equation relies on the separation of time scales: the relaxation time scale of cavity fluctuations $d_k$, on the order of $1/\kappa$, is much shorter than that of the reduced density matrix of the spins $\rho_\sigma(t)$. This separation of time scales is perfect in the steady state [@afoot; @gsc; @amitra] which we are interested in here.

The Lindblad-type operators are defined as $\mathcal{D}[X] \rho  \equiv \left( X \rho X^\dagger - X^\dagger X \rho + \rm{H.c.} \right)/2$ and the $\Gamma_{\boldsymbol{0} \to k} $’s correspond to non-equilibrium transition rates between the ground state $|\widetilde{\boldsymbol{0}} \rangle$ and a given excited many-body state $| \widetilde{k} \rangle$ \[see Eq. (\[eqqbitcollen\])\] in the one-excitation manifold of $H_\sigma$, given in Eq. (\[eq:Hspinchain\]). These are found to be $$\begin{aligned}
\label{eq:rate}
\Gamma_{\boldsymbol{0} \to k} = 2 \pi \sum_q   \Lambda_{kq}^2\rho_q(\omega_{{\rm d}}+\widetilde{E}_{\boldsymbol{0}} -\widetilde{E}_k)\,,\end{aligned}$$ where the transition matrix element is given by $$\begin{aligned}
\Lambda_{kq} =& \left| 
 \left(1+2\frac{\Delta}{\Delta_{{\rm c}}} \right) \!\!
 \frac{1}{\Delta_{{\rm q}}} \!\! \left(\frac{g}{\Delta}\right)^3 \!\!\!
 \sum_{k' k'' q'} \!\!\! f_{k k' k''} f_{q' q k'}^* \epsilon^{{\rm d}}_{k''} \epsilon^{{\rm d}}_{q'} \right|\,,\end{aligned}$$ with the tensor $f_{k k' k''} \equiv \sum_i \varphi_k(i) \varphi_{k'}^*(i) \varphi_{k''}^*(i)$. The integration over the photon-fluctuation degrees of freedom also yields Lamb-shift corrections of the energy levels, but this does not play any substantial role in our scheme, see discussion below Eq. (\[eq:wdopt\]).

#### Dynamics. {#dynamics. .unnumbered}

By virtue of having written the steady-state master equation Eq. (\[eq:master\]) in the eigenbasis of $H_{\sigma}$, all off-diagonal matrix elements of $\rho_\sigma^{\rm NESS}$ by construction vanish as the steady state is approached. Consequently, the dynamics can be faithfully described by effective rate equations for the populations of eigenstates, $n_{\boldsymbol{0}}$ and $n_k$: $$\begin{aligned}
\frac{{{\rm d}}n_{\boldsymbol{0}}}{{{\rm d}}t} &=  \gamma \, \sum_q n_q  - \Gamma_{\boldsymbol{0}\to q} \, n_{\boldsymbol{0}}  \label{eq:eom1} \\
\frac{{{\rm d}}n_k}{{{\rm d}}t} &=  - \gamma \, n_k  +   \Gamma_{\boldsymbol{0}\to k} \, n_{\boldsymbol{0}}
 + \frac{2 \gamma_\phi}{N} \sum_{q} (n_q -n_k\, ).  \label{eq:eom2}\end{aligned}$$ The terms in $\gamma$ correspond to qubit decay, flipping down the pseudo spins and relaxing the energy by $\Delta_{{\rm q}}$ (in the rotating frame). The terms in $\gamma_\phi$ correspond to pure dephasing processes, the action of which is to equalize the populations of the states in the one-excitation manifold. The emergent level structure and rates are summarized in Fig. \[fig:mediated\].

We note that, while the full dynamical evolution of the qubit-EM system (viz. Eq. (\[eq:H\]) in the presence of qubit and cavity decay) is clearly non-Markovian, the proper secularization of the equations around the operation frequency $\omega_d$ allows us to describe the qubit dynamics through the relatively transparent rate equations (\[eq:eom1\] and \[eq:eom2\]).

Irrespective of the initial conditions, Eqs. (\[eq:eom1\] and \[eq:eom2\]) have a unique non-equilibrium steady-state solution and, after transient dynamics, the occupation of the state $|k\rangle$ is given by $$\begin{aligned}
n_k^{\rm NESS} &= \frac{1}{1+2\gamma_\phi/\gamma} \frac{\Gamma_{\boldsymbol{0} \to k} + (2\gamma_\phi/N\gamma) \sum_q \Gamma_{\boldsymbol{0} \to q} }{\gamma + \sum_q \Gamma_{\boldsymbol{0} \to q}} 
\,. \label{eq:nqness}\end{aligned}$$

#### Stabilization protocol. {#stabilization-protocol. .unnumbered}

Equation (\[eq:nqness\]), together with Eq. (\[eq:rate\]), transparently elucidates how to stabilize a given pure entangled state of qubits $|k\rangle$ in the steady state. The protocol requires the maximization of $\Gamma_{\boldsymbol{0}\to k}$, given in Eq. (\[eq:rate\]), to make it the largest of all rates among ($\{ \Gamma_{\boldsymbol{0}\to q} \}$, $\gamma$, $\gamma_\phi$). This is performed by optimally tuning the drive frequency $\omega_{{\rm d}}$ such that the sharply peaked photonic spectrum $\rho_q(\omega_{{\rm d}}+\widetilde{E}_{\boldsymbol{0}} -\widetilde{E}_k)$ in Eq. (\[eq:rate\]) reaches the maximum amplitude of the Lorentzian, which is on the order of $1/\kappa$. This is possible whenever there is at least one mode $q_0$ with $\Lambda_{k q_0} \neq 0$ and the optimum $\omega_{{\rm d}}$ is the solution of the energy-conservation equation $\omega_{{\rm d}}+\widetilde{E}_{\boldsymbol{0}} -\widetilde{E}_k = \omega_{q_0}$, *i.e.* $$\begin{aligned}
\omega_{{\rm d}}= &  \frac{\omega_{{\rm q}}+ \delta\omega_{{\rm q}}+\omega_{{\rm c}}}{2} - J \cos(q_0)  \nonumber \\
& \quad + \left(\frac{g}{\Delta}\right)^2 \left[ - J \cos(k) + \frac12 \sum_{q\neq k} \frac{|\epsilon^{{\rm d}}_q|^2}{\Delta_{{\rm q}}} \right]
\,.
\label{eq:wdopt}\end{aligned}$$ This energy-matching condition describes a one-photon process in the rotating frame equivalent to a two-photon process in the laboratory frame. The corresponding Raman inelastic scattering process uses the energy of the two incoming drive photons to perform the qubit transition while simultaneously dumping a photon in one of the cavity modes [@ucbpaper]. We note that when Eq. (\[eq:wdopt\]) is satisfied, the Lamb-shift correction of $\widetilde{E}_q$ vanishes.

![\[fig:gamma\] (color online) Schematics of the pumping rate $\Gamma_{0 \to k}$ as a function of the drive frequency $\omega_{{\rm d}}$ for $N=5$. (a) Driving first cavity only: any of the five photon-fluctuation modes (responsible for the five peaks) can be used for dissipative stabilization. Driving to $|k\rangle$ while avoiding populating the nearest state $|k'\rangle$ necessitates $\kappa < \Delta \widetilde{E} \equiv |\widetilde{E}_k - \widetilde{E}_{k'}| \sim  2\pi J (g/\Delta)^2 /N$. (b) Driving all cavities equally: only one of the five modes can channel the mechanism and driving to the nearest state $|k'\rangle$ is avoided if $\kappa < \Delta \omega \equiv |\omega_k - \omega_{k'}| \sim  2\pi J / N$. ](Gamma2.eps){width="8.4cm"}

#### [Scalability and limitations.]{} {#scalability-and-limitations. .unnumbered}

Equation (\[eq:nqness\]) sets an upper bound on the fidelities, $n_{k}^{\rm NESS} \leq n^{\rm max} =
(\gamma+2\gamma_\phi/N)/(\gamma+ 2 \gamma_\phi)$, which highlights the necessity of working with qubits that have a pure dephasing rate $\gamma_\phi$ much smaller than their relaxation rate $\gamma$. This upper bound is not tight but allows us to highlight the role of the pure dephasing mechanisms. The success of the protocol, [and its scalability to large N]{}, also relies, via the numerator of Eq. (\[eq:nqness\]), on the resolving power of the spectral width of the photon density of states ($\sim \kappa$) *i.e.* the precision with which the photon fluctuations can target the spin-chain state $|k\rangle$ without exciting other eigenstates close in energy. The limitations on the resolving power depend strongly on the drive spatial profile. This is illustrated here considering two extreme cases, one where only one cavity is driven ($\epsilon^{{\rm d}}_i = \epsilon^{{\rm d}}\delta_{i1}$), the other corresponding to a case where all cavities are driven with equal amplitude ($\epsilon^{{\rm d}}_i = \epsilon^{{\rm d}}\, \forall i$). In the first case \[Fig. \[fig:gamma\](a)\], the transition rates for every $0 \to k$ have multiple peaks. Thus the stabilization of $|k\rangle$ at the optimal frequency given by Eq. (\[eq:wdopt\]) while avoiding the population of the nearest state $|k'\rangle$ requires the condition $\kappa < \Delta \widetilde{E} \equiv |\widetilde{E}_k - \widetilde{E}_{k'}| \sim  2\pi J (g/\Delta)^2 /N$. On the other hand, the second case, driving each site identically yields rates that have a single peak, that for neighboring spin chain states $k$ and $k'$ are separated by the free spectral range (of the collective EM modes) of the cavity chain, $\Delta \omega \equiv |\omega_k - \omega_{k'}| \sim  2\pi J / N$. [This indicates that the protocol with uniform driving can be scaled up to a number of qubits on the order of $N_{\rm max} \lesssim 2 \pi J / \kappa$.]{}

This situation is easily understood: the rates given in Eq. (\[eq:rate\]) contain the fulfillment of an energy-conservation condition, which can, in principle, be satisfied picking any collective EM mode $q$. The set of non-zero transition matrix elements in the sum Eq. (\[eq:rate\]) can, however, be substantially narrowed down by choosing a drive amplitude profile ($\epsilon^{{\rm d}}_{i}$) that is narrow in the momentum domain. For example, $\epsilon^{{\rm d}}_{p} = \delta_{pp_0}$, collapses the sum to a single term by imposing a quasi-momentum conservation condition between $p_0$, the target spin-chain state with momentum $k$, and $q_0$, the quasi-momentum of the collective EM mode picked for stabilization (note that conservation of momentum is strictly valid in periodic systems). This is given by $k = 2p_0 - q_0$ and is consistent with the interpretation of stabilization via a two-photon process. For such driving, the optimal frequency of the drive $\omega_{{\rm d}}$ is then given by Eq. (\[eq:wdopt\]). We note that this transparent criterion was used in Ref. [@schwartz_toward_2015] to selectively stabilize either the triplet or the singlet state of two transmon qubits.

#### Effective master equation simulations. {#effective-master-equation-simulations. .unnumbered}

To complement the analytic approach, we have performed full numerical simulations of the effective master equation (\[eq:master\]) where we (i) compute exactly the full spectrum of the spin chain $H_\sigma$ in Eq. (\[eq:Hspinchain\]) with $N=5$ qubits and open boundary conditions (ii) determine the rates between all the eigenstates and (iii) solve for the steady-state populations. In these calculations, sufficient number of higher-excitation manifolds of the spin chain were included to achieve convergence.

Because the parameter space is fairly large, we performed our simulations for a presently existing fabricated system for $N=2$ [@schwartz_toward_2015]. These parameters are quoted in the caption of Fig. \[fig:firstsiteopen\], where the fidelities to achieve various spin-chain states ($k$) for $N=5$ are compared for a localized drive \[Fig. \[fig:firstsiteopen\](a)\] and a spatially uniform drive \[Fig. \[fig:firstsiteopen\](b)\]. We note that compared to the current state of the art [@shankar_autonomously_2013], these are remarkable fidelities. These fidelities can be significantly improved by reducing the ratio of the dephasing over the qubit relaxation rate.

We have also tested the robustness of our protocol against site-to-site inhomogeneities of the different parameters and found no qualitative difference for $\delta\omega_{{\rm c}}/\omega_{{\rm c}}\sim 10^{-2}$, $\delta\omega_{{\rm q}}/\omega_{{\rm q}}\sim 10^{-4}$, $\delta g/g \sim 10^{-4}$ and $\delta J/J \sim 10^{-2}$. A more extensive analysis of achievable fidelities in the presence of site-to-site inhomogeneities will be presented in future work.

#### Periodic boundary conditions. {#periodic-boundary-conditions. .unnumbered}

Let us now consider the case of periodic boundary conditions for which the [*undriven*]{} system is space-translational and space-reversal invariant. Such a symmetry results in degeneracies between the eigenstates $|k\rangle$ and $|2\pi -k \rangle$ of the undriven spin chain (except for $k=0$ and $k=\pi$). A symmetry-breaking drive profile will generically lift the degeneracy of the spectrum and, importantly, the emergent eigenstates will strongly depend on the particular drive profile. To exemplify this point, let us start by driving the first cavity only: $\epsilon^{{\rm d}}_i = \epsilon^{{\rm d}}\delta_{i,0}$. Second-order degenerate perturbation theory in $({g}/{\Delta})(\epsilon^{{\rm d}}/\sqrt{N}\Delta_{{\rm q}})$ lifts the degeneracy in the subspaces spanned by $|k\rangle$ and $|2\pi - k \rangle$. To lowest order, the eigenstates are $$\begin{aligned}
 \label{eq:kpm}
|k_\pm \rangle \equiv \frac{| k \rangle \pm  | 2\pi - k \rangle}{\sqrt{2}}\end{aligned}$$ for all $k \in\, ]0,\pi[$ complemented with the W state $|0_+\rangle \equiv | 0 \rangle$ and $|\pi_+\rangle \equiv | \pi \rangle$ (for $N$ even), one obtains the rates $\Gamma_{\boldsymbol{0} \to k_\pm} = 2 \pi   \Lambda_{k_\pm}^2 \sum_q \rho_{q}(\omega_{{\rm d}}+\widetilde{E}_{\boldsymbol{0}} -\widetilde{E}_{k_\pm})$ with $$\begin{aligned}
 \Lambda_{k_-}  \!= 0 \mbox{ and } \Lambda_{k_+} \! =  \!\left| \frac{\sqrt{2}}{N^2}  \left(1+2\frac{\Delta}{\Delta_{{\rm c}}} \right) 
\left(\frac{g}{\Delta}\right)^3  \frac{{(\epsilon^{{\rm d}}})^2}{\Delta_{{\rm q}}} \right|
 \end{aligned}$$ for all $|k_-\rangle$ and $|k_+\rangle$ except for $|0\rangle$ or $|\pi\rangle$ in which case $\Lambda_{k_+}$ is reduced by a factor $\sqrt{2}$. Figure \[fig:firstsite\] shows that the $|k_+\rangle$ states can be obtained with substantial fidelities.

It is worth noting that, in the case of a generic driving profile $\epsilon^{{\rm d}}_i$, instead of Eq. (\[eq:kpm\]), the emergent eigenstates are given by $$\begin{aligned}
|k_\pm \rangle \equiv \frac{| k \rangle \pm  \alpha_k | 2\pi - k \rangle}{\sqrt{1+\alpha_k^2}}\end{aligned}$$ where $\alpha_k$, the relative weight of $|k\rangle$ and $|2\pi-k\rangle$, is now controlled by the ratio $\epsilon^{{\rm d}}_k/\epsilon^{{\rm d}}_{2\pi-k}$. Therefore, such a non-equilibrium symmetry-breaking scenario offers highly flexible control over the target entangled state by simply engineering the drive profile $\epsilon^{{\rm d}}_k$.

Photon-mediated interactions: General formulation
=================================================

![\[fig:3dEM\] (color online) Light-mediated interactions offer a highly versatile platform to design and control networks of interacting qubits. (a) The qubits (blue arrows) are embedded in an electromagnetic environment, which is the solution of the Maxwell equations in a given scattering geometry. (b) Integrating out the EM degrees of freedom yields effective interactions between the qubits, forming a network that can sustain large-scale entangled many-body states.](3dEM2.eps){width="7.5cm"}

\[sec:ph-mediated\] In the previous sections, we focused on a particular geometry of the cavity-qubit system, namely a one-dimensional tight-binding lattice of photons. In this section, we show that our dissipative stabilization scheme via photon-mediated interactions is broadly applicable to any engineered EM environment. The situation we consider here is depicted in Fig. 6 where the two-level systems are now placed at arbitrary locations and interact with a general EM environment \[Fig. 6(a)\]. In practice, the latter is described by the solution of Maxwell equations in an arbitrary scattering geometry, characterized by a spectral problem with certain continuity and boundary conditions. In the simplest case, this would be a single resonator or a waveguide or, as in the specific example discussed previously, an array of evanescently (or capacitively) coupled cavities (transmission line cavities). After the derivation of the most general result, we show how the effective tight-binding result can be derived from first principles.

The qubits and their coupling to the EM fluctuations are described by the Hamiltonian $$\begin{aligned}
H =  H_{\sigma}  + H_{\sigma-{\rm EM}} + H_{\rm EM}\,,\end{aligned}$$ with ($\hbar = 1$) $$\begin{aligned}
& H_{\sigma} = \sum_i \omega^{{\rm q}}_i \frac{\sigma^z_i}{2}, \quad 
H_{\rm EM} = \frac{\epsilon_0}{2} \int_\mathcal{V} \!\! {{\rm d}}^3 \boldsymbol{x} \, \left( \boldsymbol{E}^2 + c^2 \boldsymbol{B}^2 \right), \\
& H_{\sigma-{\rm EM}} =  - \int_\mathcal{V} \!\! {{\rm d}}^3 \boldsymbol{x} \,\boldsymbol{P} \cdot \boldsymbol{E}\,.\end{aligned}$$ The integrals above run over the entire volume $\mathcal{V}$ of the scattering structure. We note that such a Hamiltonian can indeed be obtained for sub-gap electrodynamics in a general superconducting circuit architecture by a proper choice of normal modes [@malekakhlagh_origin_2015], starting from the parameters (position-dependent capacitances and inductances per unit length and the parameters of the qubits) of the underlying electrical circuit.

For qubits residing at $\boldsymbol{x}_i$, each with a dipole moment strength $\mu$ (projected along a particular eigenpolarization of the electromagnetic medium), the collective atomic polarization operator can be written as $$\begin{aligned}
P(\boldsymbol{x}) = \mu \sum_i  \delta^3(\boldsymbol{x}-\boldsymbol{x}_i) \, \sigma_i^x\,.\end{aligned}$$ The electric field can generally be written in terms of a complete set of modes and corresponding eigenfrequencies $\{ \varphi_n , \omega_n \}$ specific to the chosen architecture $$\begin{aligned}
E(\boldsymbol{x}) = \sum_{n}  \varepsilon_n \varphi_n(\boldsymbol{x}) a_{n} + \mbox{H.c.}\;,\end{aligned}$$ so that $$\begin{aligned}
\label{eq:Ha1}
H_{\rm EM} = \sum_{n} \omega_n a_{n}^\dagger  a_{n}\,.\end{aligned}$$ Here, $a^\dagger_n$ ($a_n$) creates (annihilates) a photon in the spatial mode $\varphi_n(\boldsymbol{x})$ with frequency $\omega_n$ and corresponding zero-point electric field $\varepsilon_n$. The modes are assumed to satisfy the completeness and orthogonality conditions $\sum_n \varphi_n(\boldsymbol{x}) \varphi_n ^* (\boldsymbol{x}') = \delta^3(\boldsymbol{x}-\boldsymbol{x}')$ and $ \int_\mathcal{V} \! {{\rm d}}^3 \boldsymbol{x} \, \varphi_n(\boldsymbol{x}) \varphi_m^*(\boldsymbol{x}) = \delta_{nm} $. With these normalization conditions, the zero-point fields are given by $\varepsilon_n \equiv \sqrt{{\omega_n}/{2\epsilon_0}}$. Neglecting counter-rotating terms, the light-matter coupling becomes $$\begin{aligned}
H_{\sigma-{\rm EM}} = &  \sum_{i,n} g_n \varphi_n(\boldsymbol{x}_i)\,  a_n^\dagger \sigma_i^- + \mbox{H.c.}\,,\end{aligned}$$ where $\sigma_i^\pm \equiv (\sigma^x_i \pm {{\rm i}}\sigma_i^y)/2$ and $g_n \equiv - \mu \varepsilon_n$.

We shall consider the regime in which the qubit frequencies $\omega^{{\rm q}}_i$ are far detuned from the photonic modes $\omega_n$ such that the light-matter coupling can be treated via second-order perturbation theory in $g_n/(\omega^{{\rm q}}_i - \omega_n)$. This is achieved by means of a Schrieffer-Wolf (SW) transformation [@schrieffer_relation_1966] which maps $H \mapsto {{\rm e}}^{X} H  {{\rm e}}^{X^\dagger}$ where $$\begin{aligned}
 \label{eq:SWX}
X \equiv  \sum_{n,i}  \left[ \frac{g_n \varphi_n(\boldsymbol{x}_i)}{\omega_i^{{\rm q}}- \omega_n} \sigma_i^+ a_n 
-\mathrm{H.c.}
 \right]. \end{aligned}$$ This yields the following Hamiltonian $H = H_\sigma + H_{\sigma a} + H_a$ to $\mathcal{O}(g^2/\Delta^2)$, with $$\begin{aligned}
H_{\sigma} =& \sum_i \omega^{{\rm q}}_i \frac{\sigma^z_i}{2} + \frac{1}{2} \sum_{ij} \Sigma_{ij}(\omega^{{\rm q}}_i)  \sigma_j^- \sigma_i^+ + \mbox{H.c.}\,, \label{eq:hsigma} \\
H_{\sigma a} =& \sum_{i, mn}  \lambda_{i, mn}(\omega_i^{{\rm q}})  \,\frac{\sigma^z_i}{2} \, a_n a^\dagger_m   + \mbox{H.c.} \,, \label{eq:sigmaa}\end{aligned}$$ and $H_a$ is still given by Eq. (\[eq:Ha1\]). In this low-energy Hamiltonian, the second term in $H_\sigma$ describes the qubit-qubit interactions mediated by virtual photons. These photons can be emitted into and absorbed from photonic channels at frequencies $\omega_n$ with spatial distribution $\varphi_n(\boldsymbol{x})$. This is precisely the story told by the coefficients $\Sigma_{ij}(\omega) =  \sum_n |g_n|^2 \varphi_n(\boldsymbol{x}_i) \varphi_n^* (\boldsymbol{x}_j)/(\omega - \omega_n)$, which appears as a self-energy correction to the qubit sector. $\Sigma_{ij}(\omega) \sigma_j^x$ can be seen as the electric field generated at $\boldsymbol{x}_i$ by a dipole at $\boldsymbol{x}_j$, oscillating harmonically at frequency $\omega$. Note that the bare electromagnetic retarded Green’s function is given by $G^{\rm R} (\boldsymbol{x},\boldsymbol{x}'; \omega) = \sum_n \varphi_n(\boldsymbol{x}) \varphi_n^*(\boldsymbol{x}')/(\omega - \omega_n)$. This immediately implies that, in principle, all qubits interact with each other, to the extent that they can radiate EM radiation to each other \[see schematic in Fig. 6(b)\]. We note that realization-specific and restricted versions of this interaction vertex have been derived before [@majer_coupling_2007; @loo_photon-mediated_2013].

For what is proposed here, however, an equally important role is played by the term $H_{\sigma a}$ in Eq. (\[eq:sigmaa\]). This is the generalized version of the ac Stark-shift contribution to a qubit’s frequency that is well known in the dispersive regime of single-mode Cavity QED [@boissonneault_dispersive_2009], which can also be interpreted as a scattering term for photons generated by the interaction of the radiation field with qubits. This term expresses the spatial fluctuations of the effective index of refraction of the electromagnetic medium through the dynamically generated polarization fluctuations \[*i.e.* $P(\boldsymbol{x})$\] of the qubits. The interaction vertex here is again given by the resonant modes and their frequencies: $\lambda_{i,mn}(\omega) = g_n g_m^* \varphi_n (\boldsymbol{x}_i) \varphi_m^*(\boldsymbol{x}_i)/(\omega - \omega_n)$. We note that for a 1D tight-binding model of a cavity array with nearest-neighbor hopping, $\Sigma_{ij} (\omega_{{\rm q}}) \simeq  (g/\Delta)^2 [ \Delta \delta_{i,j} - J\delta_{i,j\pm 1}] $ and $\lambda_{i,kq} (\omega_{{\rm q}}) \simeq \Delta (g/\Delta)^2 \varphi_k^*(i)  \varphi_q(i)$ to lowest relevant order in $J/\Delta$ and $g/\Delta$. This result agrees with our direct derivation in Sec. \[sec:cooling\]. These results can be extended to a full-fledged stabilization protocol for a generalized W state $| W_n \rangle =  \sum_{i=0}^{N-1}  \varphi_n^* (\boldsymbol{x}_i) |\downarrow_0 \ldots\downarrow_{i-1} \, \uparrow_i \, \downarrow_{i+1} \ldots \downarrow_{N-1}  \rangle$.

Discussion and Conclusions {#sec:end}
==========================

We proposed a general and scalable method based on photon-mediated interactions to drive a set of $N$ qubits to a desired generalized W state. The particular protocol discussed here for qubits embedded in a cavity array, amounts to the dissipative stabilization of a particular [*excited state*]{} of a many-body system (in the present case, a non-integrable variant of the XY model). This approach stands in contrast to cooling techniques employed for condensed-matter and cold atomic systems that, at least in principle, target the stabilization of the ground state.

An interesting feature of the present approach compared to earlier approaches to dissipative engineering of entanglement [@kraus_preparation_2008; @diehl_quantum_2008; @cormick_dissipative_2013; @lee_emergence_2013; @rao_deterministic_2014; @reiter_scalable_2015] is the fact that [*both*]{} the unitary and the dissipative parts of the dynamics are adjusted through coupling to a common photonic bath. The Hamiltonian part provides the set of pure many-body states that can be reached in the steady state, while the dissipative part determines the occupation of those states. For cavities that are high-Q, the transitions can be made very selective. In the case of the dissipative stabilization of a W state of $N$ qubits, we discussed the scaling of the fidelity with the system size $N$.

The fact that various properties of multi-qubit dynamics can be precisely adjusted by drive parameters provides a suitable platform for quantum simulation [@rotondo_dicke_2015] and computation [@verstraete_quantum_2009]. In scaling up CQED-based simulators [[@PhysRevLett.115.240501]]{} to larger architectures, one of the main obstacles is the uncontrolled site-to-site fluctuation of system parameters [@underwood_low-disorder_2012]. In the presented scheme, the [*dynamical tuning*]{} of effective spin-chain parameters, in fact both the unitary and the dissipative parameters, through the drive frequency and amplitude provides a promising route to realize large-scale quantum simulators.

More generally, the method proposed here and its possible generalization to higher dimensional lattices holds promise for various quantum information applications, such as deterministic teleportation [@agrawal_perfect_2006; @wang_simple_2009]. The reduction of the collective dephasing mechanisms and the generation of entangled states in higher-excitation manifolds are important goals. Another interesting open question is the adaptation and extension of our protocol for targeted many-body state preparation in the photonic sector, scaling up recent approaches [@leghtas_confining_2015].

Acknowledgements {#sec:ack}
================

We are grateful to Mollie Schwartz, Leigh Martin, Emmanuel Flurin, Irfan Siddiqi, Liang Jiang, and Alexandre Blais for helpful discussions. This work has been supported by ARO Grant No. W911NF-15-1-0299 and NSF Grant No. DMR-1151810. M.K gratefully acknowledges support from the Professional Staff Congress of the City University of New York award No. 68193-0046.
