# Dissipative stabilization of entangled qubit pairs in quantum arrays with a
single localized dissipative channel

Jacopo Angeletti1,2,3, Stefano Zippilli1, David Vitali1,3,4 1 School of
Science and Technology, Physics Division, University of Camerino, I-62032
Camerino (MC), Italy
2 Department of Physics, University of Naples “Federico II”, I-80126 Napoli,
Italy
3 INFN, Sezione di Perugia, I-06123 Perugia, Italy
4 CNR-INO, I-50125 Firenze, Italy

###### Abstract

We study the dissipative stabilization of entangled states in arrays of
quantum systems. Specifically, we are interested in the states of qubits
(spin-1/2) which may or may not interact with one or more cavities (bosonic
modes). In all cases only one element, either a cavity or a qubit, is lossy
and irreversibly coupled to a reservoir. When the lossy element is a cavity,
we consider a squeezed reservoir and only interactions which conserve the
number of cavity excitations. Instead, when the lossy element is a qubit, pure
decay and a properly selected structure of $XY$-interactions are taken into
account. We show that in all cases, in the steady state, many pairs of
distant, non-directly interacting qubits, which cover the whole array, can get
entangled in a stationary way, by means of the interplay of dissipation and
local interactions.

  *

##  1 Introduction

A central, necessary ingredient of quantum technologies such as quantum
computation, simulation, and communication is the ability to control and
distribute entangled resources over large arrays of quantum systems. An
attractive strategy makes use of controlled dissipative processes to steer and
protect arrays of quantum systems into entangled states [1, 2, 3, 4, 5, 6, 7,
8, 9, 10]. In particular, it has been shown that, in order to drive the whole
system into non-trivial and potentially useful multipartite entangled states,
it is sufficient to control the dissipative dynamics of one or two localized
elements in a quantum array [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23]. This has been proved for arrays of both bosonic (by using a single
localized squeezed reservoir) [16, 17, 18, 19, 20, 21, 22] and fermionic (via
a correlated reservoir for two fermions) [15] modes. It has also been shown
that many entangled pairs can be realized by modulating the coupling between a
central cavity and two qubits (spin-1/2) in a qubit chain [23] and by
designing a correlated reservoir of two elements in an array of qubits [11,
14, 15] and cavities (bosonic modes) [11, 12, 13].

Here we present that, also in the case of qubits, it is sufficient to control
the local environment of a single element of an array to generate, in the
steady state, many entangled qubits pairs. We demonstrate this both for arrays
of cavities and qubits, and for arrays of only qubits. In the first case a
single cavity is coupled to a squeezed reservoir (see Ref. [25] for an
experimental implementation in the microwave domain and Refs. [26, 27] for its
use in the optical domain to improve the performance of gravitational wave
interferometers) and all the interactions conserve the number of excitations.
In the second, instead, a single qubit can decay, and the qubits are coupled
according to a specific geometry of $XY$-interactions. Ideally, we assume that
only one element of the arrays is lossy. We also analyze how these dynamics
are sensitive to additional noise affecting the qubits.

We note that similar states have also been found in the ground states of
certain spin Hamiltonians [28, 29, 30, 31], i.e., the so-called concentric
singlet phase [28] and rainbow states [29, 15]. They can also be generated in
a spin chain following specific dynamics [32, 33, 34], and in this context
they have been labeled nested entangled states (or matryoshka states) [32].
Moreover, analogous states are also known as thermofield double states in the
high-energy community [35, 36]. Differently from all these examples, here we
display that these states can be found as the _unique_ pure steady state of a
dissipative dynamics.

Finally, we mention the related results reported in Refs. [37, 38] where
however the entangled steady states are not unique, so that they are achieved
only if the system is prepared in peculiar initial states.

The outline of the paper is as follows. In Sec. 2 we describe in detail the
four models of arrays involving both qubits and cavities, and describe also
the main result, that is, the possibility to generate in a robust way a
stationary state of many entangled qubit pairs. In Sec. 3 we describe how a
similar dissipative generation of entangled qubit pairs can be obtained with
effective models involving only qubits. In Sec. 4 we verify our results
through the numerical solution of the dynamics of all the models presented in
the previous sections, while Sec. 5 is for concluding remarks.

##  2 Models with qubits and cavities

In this section we analyze the models which involve both qubits and cavities.
We identify four models which differ in geometries and composition of the
arrays as specified below. In all cases there is a central cavity which is
coupled to a squeezed reservoir, and the corresponding stationary state is
pure and factorized between the state of the qubits $|\psi\rangle$, and that
of the cavity/ies $|\varphi\rangle_{c}$. In particular, without exception,
each qubit gets entangled to another one, such that they form many entangled
pairs. In details, by using a proper labeling, each qubit with label $j$ gets
entangled with the $-j$th, and the steady state of the qubits can be expressed
as

$\displaystyle|\psi\rangle=\bigotimes_{j=1}^{N}\left(\sqrt{\frac{\overline{n}+1}{2\,\overline{n}+1}}\
|-\rangle_{j}|-\rangle_{-j}-\chi_{j}\
\sqrt{\frac{\overline{n}}{2\,\overline{n}+1}}\
|+\rangle_{j}|+\rangle_{-j}\right)\ ,{}$

where the variable $j$ runs over all the entangled pairs, $N$ is the number of
pairs (such that the number of qubits is $2\,N$), $|\pm\rangle_{j}$ is the
eigenvalue of the Pauli operator $\sigma_{j}^{(z)}$, for the qubit $j$, with
eigenvalue $\pm 1$, and $\overline{n}$ is the number of excitations of the
squeezed reservoir. Moreover, $\chi_{j}$ is a phase factor which depends on
the specific model [as specified in Eq. (12)]. We finally note that, by
increasing $\overline{n}$, the state of each pair in Eq. (2) tends to a Bell,
maximally-entangled state, which is a central resource in many quantum
information protocols [39].

Figure 1: (a) Chain of cavities and qubits “$Ccq$” where each cavity but the
central one interacts with a qubit. The central cavity is locally coupled to a
squeezed bath. Regardless of the initial conditions, steady-state entangled
pairs of qubits (pairs indicated by the red and green thick lines) arise. (b)
Star geometry with cavities and qubits “$Scq$”. (c) Chain of qubits with a
central cavity “$Cq$”. (d) Star of qubits with a central cavity “$Sq$”.

###  2.1 The four models

#### Chain of cavities and qubits.

The first model is an extension of the chain of cavities studied in Ref. [16]
where, here, each cavity but the central one interacts also with a qubit [see
Fig. 1 (a)]. In the following we indicate this model with the symbol “$Ccq$”,
where the upper case “$C$” stands for chain and thus addresses the geometry,
while “$cq$” the composition of the array, i.e., of cavities and qubits.

#### Star of cavities and qubits.

Similarly, the second model is the extension of a star–like bosonic array
analogous to that discussed in Ref. [17] where, here, each of the external
modes interacts with a qubit [see Fig. 1 (b)]. For this model we use the
symbol “$Scq$” where now the upper case “$S$” indicates the star geometry.

#### Chain of qubits with a central cavity.

Then, we consider a chain of qubits, where the central element is in fact a
cavity (or, in other terms, it is made of two chains coupled on one end to a
common cavity) [see Fig. 1 (c)]. In this case we use the symbol “$Cq$” for
chain of qubits.

#### Star of qubits with a central cavity.

Finally we also study a star–like array where a central cavity is coupled to
many qubits [see Fig. 1 (d)]. For this model we use the symbol “$Sq$”.

###  2.2 The master equation

In every case, the system dynamics is described by a master equation of the
form

$\displaystyle\dot{\rho}=-\frac{{\rm
i}}{\hbar}\left[H_{\xi},\rho\right]+\frac{\kappa}{2}\,{\cal L}_{c}\,\rho\ ,$
(2)

with $\xi\in\left\\{Ccq,Scq,Cq,Sq\right\\}$. The Lindblad operator ${\cal
L}_{c}$ describes dissipation via the squeezed reservoir of the central cavity
with photonic annihilation operator $b_{0}$ (see Fig. 1), and it reads

$\displaystyle{\cal L}_{c}\ \rho=2\ \beta_{0}\ \rho\
\beta_{0}^{\dagger}-\beta_{0}^{\dagger}\ \beta_{0}\ \rho-\rho\
\beta_{0}^{\dagger}\ \beta_{0}\ ,$ (3)

where

$\displaystyle\beta_{0}=\sqrt{\,\overline{n}+1}\ b_{0}-\sqrt{\,\overline{n}}\
b_{0}^{\dagger}\ ,$ (4)

is the squeezed annihilation operator. A squeezed reservoir can be realized by
driving the system with a broadband squeezed field [12, 17, 24]. In the
microwave regime, a squeezed reservoir has been experimentally realized, as
reported in Ref. [25]. This technique has also been employed to enhance the
performance of gravitational wave interferometers in the optical domain [26,
27].

###  2.3 The Hamiltonians of the four models

In all instances, we take into account only Hamiltonian interactions which
conserve the number of excitations. Thus, for two interacting cavities with
annihilation operators $b_{j}$ and $b_{k}$, we consider interaction
Hamiltonians of the form $H_{c-c}\propto b_{j}^{\dagger}\ b_{k}+h.c.$ (where
$h.c.$ indicates the Hermitian conjugate). Moreover, the cavity-qubit
interactions are described by Jaynes-Cummings Hamiltonians $H_{c-q}\propto
b_{j}^{\dagger}\,\sigma_{j}+h.c.$, with
$\sigma_{j}=\left[\sigma_{j}^{(x)}-{\rm i}\,\sigma_{j}^{(y)}\right]/2$ the
lowering operator for the qubit $j$. Finally, the interactions between two
qubits are described by an $XX$ spin-1/2 Hamiltonian
$H_{q-q}\propto\sigma_{j}^{(x)}\,\sigma_{k}^{(x)}+\sigma_{j}^{(y)}\,\sigma_{k}^{(y)}=2\,\sigma_{j}^{\dagger}\
\sigma_{k}+h.c.$.

To be specific, the Hamiltonians corresponding to each of the four models are
given by the expressions reported in Eqs. (5)-(2.3) below, which describe the
system in a reference frame rotating at the frequency of the central cavity
$\omega_{0}/2\pi$. In particular - on the one hand - in the models “$Ccq$” and
“$Scq$” with many cavities (each interacting with a qubit), all the qubits are
resonant with the central cavity, while the other cavities are detuned by a
frequency $\Delta_{c,j}$ from $\omega_{0}$. On the other hand, in the models
“$Cq$” and “$Sq$”, composed of many qubits and a central cavity, the
transition frequency of each qubit is detuned by $\Delta_{q,j}$ from
$\omega_{0}$. Our analysis relies on the assumption that $\omega_{0}$ is the
dominant parameter in the studied systems. Consequently, the frequencies of
the cavities and qubits ($\omega_{0}+\Delta_{c,j}$ and
$\omega_{0}+\Delta_{q,j}$, respectively) are orders of magnitude larger than
the introduced coupling strengths, a common feature of quantum-optical
systems. This key assumption justifies our adoption of the local master Eq.
(2) [40].

We also note that, in order to sustain the steady state of Eq. (2), the
Hamiltonians have to fulfill certain symmetry properties: they have to be
symmetric in the interaction strengths and antisymmetric in the detunings as
discussed below. This is analogous to the chiral symmetry identified in Ref.
[15] (see also Ref. [19]) and which is at the basis of the emergence of steady
state entangled pairs (equal to the ones studied here) in a finite chain of
qubits, when the two central qubits are coupled to a correlated reservoir.

Let us now introduce the explicit formulas for the Hamiltonians.

#### Chain of cavities and qubits: “$Ccq$” [Fig. 1 (a)]

$\displaystyle H_{Ccq}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\left\\{\Delta_{c,j}\left(b_{j}^{\dagger}\
b_{j}-b_{-j}^{\dagger}\
b_{-j}\right)+\left[g_{j}\left(b_{j}^{\dagger}\,\sigma_{j}+b_{-j}^{\dagger}\,\sigma_{-j}\right)\right.\right.{}$
(5)
$\displaystyle\left.\left.+\eta_{c,j}\left(b_{j}^{\dagger}\,b_{j-1}+b_{-j}^{\dagger}\,b_{-j+1}\right)+h.c.\right]\right\\}\
,$

where $j=0$ indicates the central cavity, while the positive and negative
values of $j$ indicate the elements respectively on the right and on the left
of the central cavity, with the value of $\left|{j}\right|$ measuring the
distance from the central cavity. Here we see that cavities at the same
distance on the right and on the left have opposite detuning, while all the
interactions are symmetric.

In such a situation (and also in the following), both qubits and cavities form
many entangled pairs. The dynamics of the cavities is the same as that
discussed in Ref. [16] (which specializes Ref. [22]). Here we demonstrate how
the entanglement of these bosonic modes is transferred to the qubits,
similarly to Refs. [11, 41].

#### Star of cavities and qubits: “$Scq$” [Fig. 1 (b)]

$\displaystyle H_{Scq}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\left\\{\Delta_{c,j}\left(b_{j}^{\dagger}\
b_{j}-b_{-j}^{\dagger}\
b_{-j}\right)+\left[g_{j}\left(b_{j}^{\dagger}\,\sigma_{j}+b_{-j}^{\dagger}\,\sigma_{-j}\right)\right.\right.{}$
(6) $\displaystyle\left.\left.+\eta_{c,j}\
b_{0}^{\dagger}\left(b_{j}+b_{-j}\right)+h.c.\right]\right\\}\ .$

In this regard the index $j$ does not indicate the distance from the central
cavity, but it is still used to label the elements which are entangled in the
steady state. In other words, the elements (both cavities and qubits) with
indices $j$ and $-j$ form entangled pairs. This model shares various features
with the previous one. First, the cavities are entangled in the steady state
and follow a dynamics similar to that discussed in Ref. [22] (see also Ref.
[17]) and this induces entanglement in the qubits as in the previous case [11,
41]. Second, the cavities in each pair have opposite detuning, while their
interaction coefficients are equal.

#### Chain of qubits with a central cavity: “$Cq$” [Fig. 1 (c)]

$\displaystyle H_{Cq}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\frac{\Delta_{q,j}}{2}\left[\sigma_{j}^{(z)}-\sigma_{-j}^{(z)}\right]+\hbar\left[\vphantom{\sum_{j=2}^{N}}g_{1}\
b_{0}^{\dagger}\left(\sigma_{1}+\sigma_{-1}\right)\right.{}$ (7)
$\displaystyle\left.+\sum_{j=2}^{N}\
\eta_{q,j}\left(\sigma_{j}^{\dagger}\,\sigma_{j-1}+\sigma_{-j}^{\dagger}\,\sigma_{-j+1}\right)+h.c.\right]\
,$

where, similar to the “$Ccq$” model, the positive and negative values of $j$
indicate, respectively, the qubits on the right and on the left of the central
cavity. A qubit on the right chain has a detuning which is opposite to that of
the corresponding qubit on the left chain, while the corresponding couplings
are identical. As described earlier, pairs of qubits at the same distance on
the right and on the left of the cavity are entangled in the stationary state.

#### Star of qubits with a central cavity: “$Sq$” [Fig. 1 (d)]

$\displaystyle H_{Sq}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\left\\{\frac{\Delta_{q,j}}{2}\left[\sigma_{j}^{(z)}-\sigma_{-j}^{(z)}\right]+\left[g_{j}\,b_{0}^{\dagger}\left(\sigma_{j}+\sigma_{-j}\right)+h.c.\right]\right\\}\
.{}$

As in the “$Scq$” case, the index $j$ does not indicate the distance from the
central cavity, but it is used to label the elements which are entangled in
the steady state, that is, the qubits with indices $j$ and $-j$ form entangled
pairs.

###  2.4 The steady state

Let us now analyze in detail the steady state of the models defined above.
When $g_{j}=0$ (indicating no interaction with the qubits), the squeezed
reservoir drives the central cavity towards a squeezed state. Correspondingly,
as shown in Refs. [16, 17, 22], the other cavities (in the models “$Ccq$” and
“$Scq$”) approach a pure entangled state constituted of many two-mode squeezed
states. It can be expressed as

$\displaystyle|\,\varphi\rangle_{c}=U_{c}|0\rangle_{c},$ (9)

where $|0\rangle_{c}$ is the vacuum and $U_{c}$ the unitary which generates
the steady state. In detail, for the models of cavities and qubits (“$Ccq$”
and “$Scq$”), $U_{c}$ is the product of a squeezing operator on the central
cavity mode, and many two-mode squeezing operators for the modes with opposite
indices

$U_{c}={\rm e}^{\frac{r}{2}\left({b_{0}^{\dagger}}^{2}-b_{0}^{2}\right)}\
\bigotimes_{j=1}^{N}{\rm
e}^{\chi_{j}\,r\,\left(b_{j}^{\dagger}\,b_{-j}^{\dagger}-b_{j}\,b_{-j}\right)}\
,\hskip 14.22636pt{\rm for}\ \ Ccq\ \ {\rm and}\ \ Scq\ ,$ (10)

where $\tanh(r)=\sqrt{\,\overline{n}/\left(\,\overline{n}+1\right)}$. Instead,
for the models “$Cq$” and “$Sq$”, $U_{c}$ is the single mode squeezing
operator

$\displaystyle U_{c}={\rm
e}^{\frac{r}{2}\left({b_{0}^{\dagger}}^{2}-b_{0}^{2}\right)}\ ,\hskip
14.22636pt{\rm for}\ \ Cq\ \ {\rm and}\ \ Sq\ .$ (11)

In particular the factor $\chi_{j}$, which appears in Eqs. (2) and (10), is
given by

$\chi_{j}=\left\\{\begin{array}[]{llll}&(-1)^{j}\ ,&&{\rm for}\ \ Ccq\ ,\\\
&(-1)^{j+1}\ ,&&{\rm for}\ \ Cq\ ,\\\ &-1\ ,&&{\rm for}\ Scq\ ,\\\ &1\ ,&&{\rm
for}\ Sq\end{array}\ ,\\\ \right.$ (12)

for $j\neq 0$. We also note that these operators [Eqs. (10) and (11)] realize
the Bogoliubov transformation

$\displaystyle U_{c}^{\dagger}\ b_{j}\ U_{c}$ $\displaystyle=$
$\displaystyle\sqrt{\,\overline{n}+1}\ b_{j}+\chi_{j}\ \sqrt{\,\overline{n}}\
b_{-j}^{\dagger}\ ,$ (13)

for all $j$, also for $j=0$ with $\chi_{0}=1$.

Now, one can check that, in general, the product state

$\displaystyle|\Psi\rangle=|\,\varphi\rangle_{c}\ |\psi\rangle\ ,$ (14)

between the state of the qubits given by Eq. (2), and that of the cavity/ies
given by Eq. (9), with $U_{c}$ and $\chi_{j}$ defined in Eqs. (10)-(12), _is a
steady state for the four models_. In other words, one finds that $-\frac{{\rm
i}}{\hbar}\left[H_{\xi},\left|\Psi\right\rangle\left\langle\Psi\right|\right]+\frac{\kappa}{2}\,{\cal
L}_{c}\,\left|\Psi\right\rangle\left\langle\Psi\right|=0$.

This is the result of the destructive interference which takes place when
these systems are endowed with the specific symmetries described in Sec. 2.3.
To be specific, by dividing the Hamiltonians as the sum of the term which
involves only the cavity operators $H_{c,\xi}$ (this is non-zero only for the
models “$Ccq$” and “$Scq$”), the one for the qubits alone $H_{q,\xi}$ (this is
non-zero only for the models “$Cq$” and “$Sq$”) and the terms which describe
the interactions between cavity/ies and qubits $H_{c-q,\xi}$, such that

$\displaystyle H_{\xi}=H_{c,\xi}+H_{q,\xi}+H_{c-q,\xi}\ ,$ (15)

we find that

$\displaystyle-\frac{{\rm
i}}{\hbar}\left[H_{c,\xi},\left|\,\varphi\right\rangle_{c\\!}\left\langle\varphi\right|\right]+\frac{\kappa}{2}\,{\cal
L}_{c}\,\left|\,\varphi\right\rangle_{c\\!}\left\langle\varphi\right|=0\ ,$
(16)

as demonstrated in Refs. [16, 19, 22], and

$\displaystyle H_{q,\xi}\,|\psi\rangle$ $\displaystyle=$ $\displaystyle 0\ ,$
(17) $\displaystyle H_{c-q,\xi}\,|\Psi\rangle$ $\displaystyle=$ $\displaystyle
0\ ,$ (18)

because of the destructive quantum interference between transitions which
involve the qubit states in the quantum superposition of Eq. (2).

###  2.5 Dynamics in the squeezed representation

Figure 2:  Low-excitations eigenlevels (black horizontal lines) corresponding
to the Hamiltonian $\left.{\widetilde{H}_{Ccq}}\right|_{\eta_{c,j}=0,g_{j}=0}$
(in the squeezed representation) without interaction terms [see Eqs. (5), (21)
and (2.5)] and corresponding matrix elements of the interaction terms (red,
blue, green, and orange arrows; each color marks a different coupling strength
as shown in the right-bottom corner of the figure), for the model “$Ccq$” with
$N=1$ (which is equal to the model “$Scq$” with $N=1$). The purple wavy arrow
illustrates the transfer of population due to the decay of the central cavity;
$|\psi\rangle=\left(\sqrt{\,\overline{n}+1}\ |--\rangle+\sqrt{\,\overline{n}}\
|++\rangle\right)/\sqrt{2\,\overline{n}+1}$ is the qubit state of Eq. (2) and
$|\psi^{\perp}\rangle=\left(\sqrt{\,\overline{n}}\
|--\rangle-\sqrt{\,\overline{n}+1}\
|++\rangle\right)/\sqrt{2\,\overline{n}+1}$ is its orthogonal state. The
parameters $n_{j}$ and $n\equiv\sum_{j}\,n_{j}$ indicate the Fock states
$|n_{-1}\ n_{0}\ n_{1}\rangle_{c}$ of the cavities in the squeezed
representation (i.e., the squeezed Fock states $U_{c}|n_{-1}\ n_{0}\
n_{1}\rangle_{c}$ in the original representation). The values of $\lambda$ tag
the eigenvalues of $\left.{\widetilde{H}_{Ccq}}\right|_{\eta_{c,j}=0,g_{j}=0}$
corresponding to the various group of levels. The horizontal dashed gray lines
divide the groups of levels associated to each eigenvalue. Levels
corresponding to the same eigenvalue (in the same group) are reported at
different vertical positions in order to provide a better visualization. The
state $|\psi\rangle$ with zero excitations, $n=0$, is the steady state of the
system dynamics, where the population accumulates. In fact, being the only
state decoupled from all the other levels, it does not lose population but it
is populated by the cavity decay.

In order to verify Eq. (18) and gain insight into the steady state dynamics,
it is useful to analyze the system in the representation in which the cavity
steady state is the vacuum. Namely, we consider the master equation for the
transformed density matrix $\widetilde{\rho}=U_{c}^{\dagger}\,\rho\,U_{c}$,
which is given by

$\displaystyle\dot{\widetilde{\rho}}=-\frac{{\rm
i}}{\hbar}\left[\widetilde{H}_{\xi},\widetilde{\rho}\right]+\frac{\kappa}{2}\,\widetilde{\cal
L}_{c}\,\widetilde{\rho}\ ,$ (19)

with the Lindblad term which describes dissipation of the zeroth mode being

$\displaystyle\widetilde{\cal L}_{c}\ \widetilde{\rho}=2\ b_{0}\
\widetilde{\rho}\ b_{0}^{\dagger}-b_{0}^{\dagger}\ b_{0}\
\widetilde{\rho}-\widetilde{\rho}\ b_{0}^{\dagger}\ b_{0}\ .$ (20)

The transformed Hamiltonian in Eq. (19) can be written as

$\displaystyle\widetilde{H}_{\xi}$ $\displaystyle=$ $\displaystyle
U_{c}^{\dagger}\ H_{\xi}\ U_{c}{}$ (21) $\displaystyle=$ $\displaystyle
H_{c,\xi}+H_{q,\xi}+\widetilde{H}_{c-q,\xi}\ ,$

where the Jaynes-Cummings interaction term
$\widetilde{H}_{c-q,\xi}=U_{c}^{\dagger}\ H_{c-q,\xi}\ U_{c}$ may be expressed
as

$\displaystyle\widetilde{H}_{c-q,\xi}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\ g_{j}\left(b_{j}^{\dagger}\
\tau_{j}+b_{-j}^{\dagger}\ \tau_{-j}+h.c.\right)\ ,\hskip 14.22636pt{\rm for}\
\xi\in\left\\{Ccq,Scq\right\\}\ ,{}$ $\displaystyle\widetilde{H}_{c-q,Cq}$
$\displaystyle=$ $\displaystyle\hbar\ g_{1}\
b_{0}^{\dagger}\left(\tau_{1}+\tau_{-1}\right)+h.c.\ ,{}$
$\displaystyle\widetilde{H}_{c-q,Sq}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\
g_{j}\left[b_{0}^{\dagger}\left(\tau_{j}+\tau_{-j}\right)+h.c.\right]\ ,$ (22)

with $\tau_{j}$ given by the collective qubit operator

$\displaystyle\tau_{j}=\sqrt{\,\overline{n}+1}\ \sigma_{j}+\chi_{j}\
\sqrt{\,\overline{n}}\ \sigma_{-j}^{\dagger}\ .$ (23)

Now it is easy to verify that the transformed state

$\displaystyle\widetilde{|\Psi\rangle}=\widetilde{|\,\varphi\rangle}_{c}\
|\psi\rangle=|0\rangle_{c}\ |\psi\rangle\ ,$ (24)

(with $\widetilde{|...\rangle}=U_{c}^{\dagger}\,|...\rangle$) fulfills the
relation $\widetilde{H}_{c-q,\xi}\,\widetilde{|\Psi\rangle}=0$ [which is
equivalent to Eq. (18)]. In fact, on the one hand,
$b_{j}\,\widetilde{|\Psi\rangle}=0$ because the cavity(-ies) is (are) in the
vacuum (in this representation) and, on the other hand,

$\displaystyle\tau_{j}\,|\psi\rangle=0\ ,\hskip 28.45274pt\forall\ j\ .$ (25)

In other terms, in this representation, all the cavities dissipate (through
the central cavity) and approach the vacuum. Correspondingly the qubits
population accumulates in the entangled state of Eq. (2) in a way similar to
optical pumping. In fact, the state of Eq. (24) is the only one which does not
decay and it is decoupled from all the others, and in the meanwhile it is
populated by the decay of the central cavity.

It is instructive to visualize this with a simple example, i.e., when $N=1$,
for the models “$Ccq$” and “$Scq$” which are equal. In Fig. 2 we report the
eigenstates of the system Hamiltonian without interactions and use arrows to
connect levels which are coupled by the interaction terms. The number of
cavity excitations (in the squeezed representation) increases from left to
right, and the dissipation of the central cavity is responsible for an
irreversible transfer of population from the levels on the right to those on
the left. The figure shows that the qubit state of Eq. (24) [i.e., the state
of Eq. (2) with zero cavity excitations] is the only state which is decoupled
from the other levels, and at the same time it does not decay and is populated
by the cavity decay. As a consequence this state is stable in the steady
state. Similar considerations hold also for the other models.

##  3 Models with only qubits

A similar qubits dynamics (in the squeezed representation) is observable also
if we modify our models by replacing all the cavity modes with fresh new
qubits. Namely, we may consider the master Eq. (19) and replace all the
$b_{j}$ and $b_{j}^{\dagger}$ with other lowering and rising qubits operators
$\sigma_{c,j}$ and $\sigma_{c,j}^{\dagger}$ (where the label c indicates that
these are the new qubits in place of the cavities of Sec. 2 and Fig. 1). In
this way we find new models which consist of only qubits, with a single
central lossy one. Moreover, the qubits interact with a peculiar structure
(different for each model) of $XY$-interactions (see Fig. 3).

Figure 3:  Sketch of the arrays of qubits analogous of the models of Fig. (1)
but with the cavities replaced by additional qubits. In each model the central
qubit is lossy (described by the wavy line) and the qubits which are entangled
in the steady state are indicated by the red and green thick lines. In (a) and
(b) the solid black lines connecting two qubits indicate $XX$\- (isotropic
$XY$-) interactions
$\propto\sigma_{j}^{(x)}\sigma_{k}^{(x)}+\sigma_{j}^{(y)}\sigma_{k}^{(y)}$,
while the dotted black lines indicate interactions of the form
$\propto\sigma_{j}^{(x)}\sigma_{k}^{(x)}-\sigma_{j}^{(y)}\sigma_{k}^{(y)}$. In
(c) and (d) the interactions between the central spin and the neighboring ones
is of $XY$-type with anisotropic couplings $\propto
g^{(x)}\,\sigma_{j}^{(x)}\sigma_{k}^{(x)}+g^{(y)}\,\sigma_{j}^{(y)}\sigma_{k}^{(y)}$,
while the black lines in (c) between the other spins indicate
$XX$-interactions
$\propto\sigma_{j}^{(x)}\sigma_{k}^{(x)}+\sigma_{j}^{(y)}\sigma_{k}^{(y)}$.

Let us explicitly write down the equations for these models. The master
equation has the form of Eq. (19)

$\displaystyle\dot{\varrho}=-\frac{{\rm
i}}{\hbar}\left[H_{\xi}^{\circ},\varrho\right]+\frac{\kappa}{2}\,{\cal
L}_{q}\,\varrho\ ,$ (26)

where also in this case we use the labels
$\xi\in\left\\{Ccq,Scq,Cq,Sq\right\\}$ to highlight the relation with the
models of Fig. 1, but where now the Hamiltonian and the Lindblad operator
include only qubits operators, following the substitution
$b_{j}\to\sigma_{c,j}$. In other words

$\displaystyle{\cal L}_{q}$ $\displaystyle=$
$\displaystyle\left.{\widetilde{\cal L}_{c}}\right|_{b_{j}\to\sigma_{c,j}}\
,{}$ $\displaystyle H_{\xi}^{\circ}$ $\displaystyle=$
$\displaystyle\left.{\widetilde{H}_{\xi}}\right|_{b_{j}\to\sigma_{c,j}}\ ,$
(27)

namely

$\displaystyle{\cal L}_{q}\ \varrho=2\ \sigma_{c,0}\ \varrho\
\sigma_{c,0}^{\dagger}-\sigma_{c,0}^{\dagger}\ \sigma_{c,0}\ \varrho-\varrho\
\sigma_{c,0}^{\dagger}\ \sigma_{c,0}\ ,$ (28)

and

$\displaystyle H_{\xi}^{\circ}$ $\displaystyle=$ $\displaystyle
H_{c,\xi}^{\circ}+H_{q,\xi}+H_{c-q,\xi}^{\circ}\ ,$ (29)

where, as in Sec. 2, $H_{q,\xi}$ is non-zero only for
$\xi\in\left\\{Cq,Sq\right\\}$ and $H_{c,\xi}^{\circ}$ is non-zero only for
$\xi\in\left\\{Ccq,Scq\right\\}$, with

$\displaystyle H_{c,Ccq}^{\circ}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\sum_{\zeta=\pm}\left[\zeta\Delta_{c,j}\
\sigma_{c,\zeta j}^{(z)}+\frac{\eta_{c,j}}{2}\left(\sigma_{c,\zeta
j}^{(x)}\,\sigma_{c,\zeta(j-1)}^{(x)}+\sigma_{c,\zeta
j}^{(y)}\,\sigma_{c,\zeta(j-1)}^{(y)}\right)\right]\ ,{}$ $\displaystyle
H_{c,Scq}^{\circ}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\sum_{\zeta=\pm}\left[\zeta\Delta_{c,j}\
\sigma_{c,\zeta
j}^{(z)}+\frac{\eta_{c,j}}{2}\left(\sigma_{c,0}^{(x)}\,\sigma_{c,\zeta
j}^{(x)}+\sigma_{c,0}^{(y)}\,\sigma_{c,\zeta j}^{(y)}\right)\right]\ ,{}$

and moreover the interaction terms, derived from the Jaynes-Cummings terms of
the previous section, are given by the following expressions. Indeed, in the
models “$Ccq$” and “$Scq$”, each qubit corresponding to a cavity of the
previous model interacts with two qubits according to the Hamiltonians

$\displaystyle H_{c-q,\xi}^{\circ}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\frac{g_{j}}{2}\sum_{\zeta=\pm}\left\\{\sigma_{c,\zeta
j}^{(x)}\left[\sqrt{\,\overline{n}+1}\,\sigma_{\zeta
j}^{(x)}+\chi_{j}\,\sqrt{\,\overline{n}}\,\sigma_{-\zeta
j}^{(x)}\right]\right.{}$ $\displaystyle\left.+\sigma_{c,\zeta
j}^{(y)}\left[\sqrt{\,\overline{n}+1}\,\sigma_{\zeta
j}^{(y)}-\chi_{j}\,\sqrt{\,\overline{n}}\,\sigma_{-\zeta
j}^{(y)}\right]\right\\}\ ,\hskip 14.22636pt{\rm for}\
\xi\in\left\\{Ccq,Scq\right\\}\ .{}$

Instead, in the models “$Cq$” and “$Sq$”, the Jaynes-Cummings terms result in
the anistrotopic $XY$-interaction Hamiltonians

$\displaystyle H_{c-q,Cq}^{\circ}$ $\displaystyle=$
$\displaystyle\hbar\sum_{\zeta=\pm}\left[g_{1}^{(x)}\ \sigma_{c,0}^{(x)}\
\sigma_{\zeta 1}^{(x)}+g_{1}^{(y)}\ \sigma_{c,0}^{(y)}\ \sigma_{\zeta
1}^{(y)}\right]\ ,{}$ $\displaystyle H_{c-q,Sq}^{\circ}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j}\sum_{\zeta=\pm}\left[g_{j}^{(x)}\
\sigma_{c,0}^{(x)}\ \sigma_{\zeta j}^{(x)}+g_{j}^{(y)}\ \sigma_{c,0}^{(y)}\
\sigma_{\zeta j}^{(y)}\right]\ ,{}$

with

$\displaystyle g_{j}^{(x)}$ $\displaystyle=$
$\displaystyle\frac{g_{j}}{2}\left(\sqrt{\,\overline{n}+1}+\chi_{j}\,\sqrt{\,\overline{n}}\right)\
,{}$ $\displaystyle g_{j}^{(y)}$ $\displaystyle=$
$\displaystyle\frac{g_{j}}{2}\left(\sqrt{\,\overline{n}+1}-\chi_{j}\,\sqrt{\,\overline{n}}\right)\
,$ (30)

where $\chi_{j}$ is defined in Eq. (12).

Now it is easy to check that, as in the previous section, a steady state for
these models is

$\displaystyle|\Psi^{\circ}\rangle=|-\rangle_{c}\ |\psi\rangle\ ,$ (31)

where $|-\rangle_{c}$ indicates the state for all the qubits with lowering
operator $\sigma_{c,j}$, where each qubit is in the eigenstate of
$\sigma_{c,j}^{(z)}$ with eigenvalue $-1$.

##  4 Numerical results

The fact that the steady state Eq. (14) is in fact unique, for specific
choices of the parameters $\Delta_{c,j}$, $\eta_{c,j}$, $\Delta_{q,j}$,
$\eta_{q,j}$, and $g_{j}$, can be numerically verified. In this section we
report the numerical results evaluated in the squeezed representation for
various sizes of the arrays of the four models of Sec. 2, and by truncating
the Hilbert space of the cavities at various Fock numbers. We also include the
limiting case of only two levels for the cavity modes, which hence corresponds
to the qubits models of Sec. 3.

Additionally, we investigate the sensitivity of these dynamics to extra noise,
as depicted in Figs. 4-7, and examine their scaling behavior with the size of
the array, as shown in Fig. 7. We consider Eq. (19) and include phase noise on
the qubits 111The effect of additional cavity decay on similar systems (made
only of cavities) has been analyzed in detail in Refs. [16, 22], showing that
the entanglement dynamics is unaffected as far as the total additional decay
rate is smaller than the coupling rate to the squeezed reservoir $\kappa$.,
according to the equation $\dot{\widetilde{\rho}}{}={\cal L}\
\widetilde{\rho}$, where the total Liouvillian superoperator is given by

$\displaystyle{\cal L}\ \widetilde{\rho}=-\frac{{\rm
i}}{\hbar}\left[\widetilde{H}_{\xi},\widetilde{\rho}\right]+\frac{\kappa}{2}\,\widetilde{\cal
L}_{c}\,\widetilde{\rho}+\gamma\,{\cal D}\,\widetilde{\rho}\ ,$ (32)

and where the additional noise on the qubits is described by the Lindblad term

$\displaystyle{\cal D}\
\widetilde{\rho}=\sum_{\mbox{\scriptsize$\begin{array}[]{c}j=-N\\\ j\neq
0\end{array}$}}^{N}\ \sigma_{j}^{(z)}\ \widetilde{\rho}\
\sigma_{j}^{(z)}-\widetilde{\rho}\ .$ (35)

We highlight that the structure of the Hamiltonians discussed in the previous
sections is not sufficient to guarantee that the steady state is unique. A
trivial example is when, in the model “$Sq$”, all the detunings and all the
couplings are equal, such that $\Delta_{q,j}=\Delta_{q,j^{\prime}}$ and
$g_{j}=g_{j^{\prime}}$ for all $j,j^{\prime}\in\left\\{1,\dots,\,N\right\\}$.
In this highly symmetric case, any partition of the qubits in pairs may be
used to construct a state as that of Eq. (2) which is stationary. All these
possible partitions will give rise to a stationary subspace of the system
dynamics. However, the actual steady state will depend on the initial state
and typically be a statistical mixture of states in this subspace. Therefore,
to obtain a unique steady state, such as the pure steady state discussed in
Sec. 2.4 (or an approximate mixed-state version thereof in the presence of
finite $\gamma$), it is necessary to avoid these highly symmetric situations,
e.g, by using different values of detunings and couplings for each pair.

Here, we characterize the steady state in terms of the concurrence [39]
between pairs of qubits. In particular, we verify numerically that, for the
chosen set of parameters, only qubits with opposite indices $j$ and $-j$ can
get entangled in the steady state.

Figure 4: Time evolution of the concurrence for the models: $``Cq"$ with a
single qubit pair, $N=1$, in (a) and (d), $``Cq"$ with $N=2$ in (b) and (e),
and $``Ccq"$ with $N=1$ in (c) and (f), and for various values of the
dimension of the working Hilbert space of the cavities (in the squeezed
representation) as reported in the legend. In all the plots $\overline{n}=1$.
In (b) and (e) the solid (dashed) lines refer to the couple $j=2\,(1)$. In the
first row [(a), (b), and (c)] the dephasing rate is $\gamma=0$. In the second
row [(d), (e), and (f)] $\gamma=5\times 10^{-4}\kappa$. In (a) and (d)
$\Delta_{q,1}\simeq-0.193\kappa$ and $g_{1}=0.36\kappa$. In (b) and (e)
$\Delta_{q,1}\simeq-0.193\kappa$, $\Delta_{q,2}=\Delta_{q,1}+0.05\kappa$,
$g_{1}=0.36\kappa$, and $\eta_{q,1}=0.362\kappa$. In (c) and (f), the
dimensions of the Hilbert spaces of all the cavities are truncated at the same
Fock number $n_{0}^{(\rm max)}=n_{\pm 1}^{(\rm max)}$, and
$\Delta_{c,1}\simeq-0.26\kappa$, $g_{1}=0.36\kappa$. The values of the
detunings are chosen in order to maximize the decay rate of the arrays for the
given choice of interaction strengths (which here are chosen for simplicity of
the same order of magnitude). This is done by maximizing the real part of the
smallest non-zero value of the total Liouvillian ${\cal L}$ as a function of
$\Delta$, see Fig. 5 (c) and Figs. 6 (c), (f), and (d) for specific examples.

We compute the evolution of the system numerically with wave function Monte
Carlo techniques, and analyze the stability of the result by truncating the
Hilbert space of the cavities to various Fock numbers $n_{j}^{(\rm max)}$ (see
Fig. 4). To ensure the uniqueness of the steady state of the simulated models,
we perform numerical analysis on the spectrum of the total Liouvillian. Our
results confirm that, for the chosen parameters, the steady state is indeed
unique, as evidenced by the null space of the Liouvillian having dimension
one. Due to the numerical complexity of the problem, we only consider small
arrays. In the squeezed representation one can use a relatively low number of
Fock states, with the lowest value $n_{j}^{(\rm max)}=1$ corresponding to
models with only qubits of Sec. 3 and the largest values of $n_{j}^{(\rm
max)}$ which approach the models which include the cavities of Sec. 2. In the
case of qubit-only models, we can simulate the full Hilbert space, which
allows us to unambiguously demonstrate the uniqueness of the steady state for
the chosen parameters. However, for models involving cavities, one can only
simulate a finite part of the infinite-dimensional Hilbert space, raising
concerns about the effective uniqueness of the steady state for the complete
models. Nevertheless, an analysis of Fig. 2 suggests that there are no
subspaces outside of the simulated part that do not dissipate and are
disconnected from it. This indicates that the steady state is likely to be
unique for these cases as well.

Figure 5: (a) Steady state concurrence for the model $``Cq"$ with $N=1$ (that
is equal to $``Sq"$ with $N=1$) as a function of the average number of
excitations of the squeezed reservoir $\overline{n}$, and for various values
of the dephasing rate $\gamma$. (b) Time evolution of the concurrence for the
highlighted (yellow) points in (a) corresponding to $\overline{n}=0.5$ (dark
red), $1$ (red), and $5$ (salmon) with $\gamma=5\times 10^{-4}\kappa$. (c)
Real part of the first two eigenvalues of $\mathcal{L}$ for the same points.
The initial state is vacuum (in the squeezed representation) for the cavity
and the eigenstate of $\sigma_{j}^{(z)}$ with eigenvalue $-1$ for all the
qubits. The vertical black, dashed line in (c) indicates the value of
$\Delta_{q,1}\simeq 0.197\kappa$. The coupling strength is $g_{1}=0.36\kappa$.
The cavity Fock space is truncated at $n_{0}^{(\rm max)}=2$ (in the squeezed
representation) to approximate a hybrid model of cavities and qubits. We have
verified that similar results hold also for models of only qubits
($n_{0}^{(\rm max)}=1$) and for larger values of $n_{0}^{(\rm max)}$. (d)
Decay rate analysis: same results of (b) where we subtract - at each curve -
its steady state value, take the modulus, and fit the result with an
exponential decay. The slope of these lines indicate the rate of decay towards
the steady state. These rates are consistent withe the real part of the
eigenvalues of ${\cal L}$ identified by the vertical dashed line in plot (c).

In Fig. 4 (a)-(c) we observe that in the ideal case ($\gamma=0$) the steady
state concurrence is the same for all the models, and independent from the
dimension of the Hilbert space. This confirms the result of Eqs. (2) and (14),
i.e., that without any additional dissipation channel the steady state of the
qubits is the same for all the models and it depends only on the amount of
squeezing of the reservoir which is determined by $\bar{n}$. However, we
notice that the dynamics involving cavities [corresponding to larger values of
$n_{j}^{(\rm max)}$] is significantly faster than that of the models with
qubits only [$n_{j}^{(\rm max)}=1$]. On the contrary, dephasing reduces the
final entanglement, and maximum reduction is observed for the slowest models,
namely the models with only qubits [see Figs. 4 (d)-(f)].

Figure 6: Results similar to those of Fig. (5) (a)-(c) for the models “$Cq$”
with $N=2$ (a)-(c), “$Sq$” with $N=2$ (d)-(f), and “$Ccq$” with $N=1$ (that is
equal to “$Scq$” with $N=1$) (g)-(i). In (a) and (d) the solid (dashed) lines
refer to the couple $j=2\,(1)$. In (a)-(c) $\Delta_{q,1}\simeq-0.13\kappa$,
$\Delta_{q,2}=\Delta_{q,1}+0.05\kappa$, $g_{1}=0.36\kappa$, and
$\eta_{q,1}=0.362\kappa$. The Fock space of the cavity is truncated at
$n_{0}^{(\rm max)}=2$. In (d)-(f) $\Delta_{q,1}\simeq 1.3\kappa$,
$\Delta_{q,2}=\Delta_{q,1}+0.05\kappa$, $g_{1}=0.36\kappa$, and
$g_{2}=0.362\kappa$. The Fock space of the cavity is truncated at $n_{0}^{(\rm
max)}=2$. In (g)-(i) $g_{1}=0.36\kappa$ and $\Delta_{q,1}\simeq-0.17\kappa$.
The Fock spaces of the cavities are truncated at $n_{0}^{(\rm max)}=2$ and
$n_{\pm 1}^{(\rm max)}=1$.  Figure 7: (a) Time evolution of the concurrence of
the $j$th couple for different numbers of qubit pairs $N\leq 5$, for “$Cq$”.
$N=1$ purple, $N=2$ blue, $N=3$ green, $N=4$ orange, and $N=5$ red.
$\overline{n}=1$ and $\gamma/\kappa=10^{-5}$. The different qubit pairs are
denoted with different symbols. Cross: $|j|=1$ pair; triangle: $|j|=2$ pair;
circle: $|j|=3$ pair; diamond: $|j|=4$ pair; square: $|j|=5$ pair. Only the
results for the first pair is suitably coloured, while all the others are
gray-scale. The Fock space of the central cavity is truncated to $n_{0}^{(\rm
max)}=2$. The system parameters are $g_{1}=0.36\kappa$,
$\eta_{q,j}=g_{1}+0.002\kappa(j-1)$, and
$\Delta_{q,j}=\Delta^{\star}+0.05\kappa(j-1)$, where the value of
$\Delta^{\star}$ is different for each $N$ as reported in the legend, and it
has been chosen in order to maximize the concurrence. (b) Corresponding steady
state concurrence for all the pairs.

Figs. 5 and 6 show that, on the one hand, at zero additional noise
($\gamma=0$) maximum entanglement is achieved for larger values of
$\overline{n}$, as expected from Eq. (2), which approaches the product of many
Bell, maximally entangled states. On the other hand, when the rate $\gamma$ of
the additional noise is finite, maximum entanglement is achieved at _finite_
values of $\overline{n}$. In fact, the larger $\overline{n}$, the slower the
dynamics, as illustrated by these figures [see the decay rates for different
values of $\bar{n}$ in Fig. 5 (d)] and, as a consequence, if $\gamma$ is too
large noise and decoherence have enough time to spoil the generation of the
steady state.

The fact that a larger $\bar{n}$ corresponds to a slower dynamics is described
by Figs. 5 (b), (c), and (d) (and the corresponding plots in Fig. 6).
Specifically, in Fig. 5 (c) [and Figs. 6 (c), (f), and (i)] we report the real
part of the eigenvalues of the total Liouvillian ${\cal L}$ [see Eq. (32)].
The smallest (in modulus) real part determines the rate of decay towards the
steady state: a larger (in modulus) real part corresponds to a faster
dynamics. This is clearly described by Fig. 5 (d), where we report the time
evolution of the concurrence relative to its steady state value. This quantity
describes an exponential decay with rate given by the values of the
eigenvalues identified in Fig. 5 (c). We verified this behavior also for the
other models.

In Fig. 7 we verify numerically that the steady state remains unique even for
larger number of qubit pairs, when the values of the detuning and of the
couplings are properly selected. We show the concurrence for the model “$Cq$”
with up to $10$ qubits. Both the detuning and the couplings vary linearly with
the pair index according, respectively, to the relations
$\Delta_{q,j}=\Delta^{\star}+0.05\kappa(j-1)$ and
$\eta_{q,j}=g_{1}+0.002\kappa(j-1)$.

In this figure, we also maximize the final concurrence as a function of
$\Delta^{\star}$ only. Specifically, the value of $\Delta^{\star}$ is chosen
in order to maximize the smallest real part of the eigenvalues of ${\cal L}$
for each $N$ (when all the other parameters are kept fixed), which – as
discussed above – determines how fast the system approaches the steady state.
These results show that the final concurrence decreases with the number of
pairs and display a behaviour very similar to the entanglement obtained in
chains of bosonic modes [16].

Figure 8:  Steady state concurrence for the model “$Cq$” of only qubits with
$N=1,\,2$, and $3$ entangled pairs, evaluated for random variations of the
system parameters from the symmetric configuration used in the previous
figures. The results in the left column are evaluated relaxing the condition
of opposite detunings for qubits with opposite indices: in place of
$\Delta_{q,j}$ we use $\Delta_{q,\pm j}^{(r)}=\pm\Delta_{q,j}\left(1+r_{\pm
j}^{(\Delta)}\right)$ with $r_{\pm j}^{(\Delta)}$ random variables uniformly
distributed in the range $\left[-r_{\rm max},r_{\rm max}\right]$. The results
in the right column are evaluated relaxing the condition of equal couplings
for qubits with opposite indices: in place of $g_{1}$ and $\eta_{j}$ we use
$g_{\pm 1}^{(r)}=g_{1}\left(1+r_{\pm j}^{(g)}\right)$ and $\eta_{\pm
j}^{(r)}=\eta_{j}\left(1+r_{\pm j}^{(\eta)}\right)$, respectively, with
$r_{\pm 1}^{(g)}$ and $r_{\pm j}^{(\eta)}$ random variables uniformly
distributed in the range $\left[-r_{\rm max},r_{\rm max}\right]$. In each
plot, for each value of $r_{\rm max}$, we report $200$ random realizations
(gray lines), the corresponding averages (red dots), and standard deviations
(red lines). The blue lines are the ideal results evaluated for the symmetric
configuration with $r_{j}^{(\Delta)}=r_{j}^{(g)}=r_{j}^{(\eta)}=0$. In all
plots $\Delta_{j}=\left(1.8+0.2\,j\right)\kappa$, $g_{1}=0.3\,\kappa$,
$\eta_{j}=\left(0.45+0.05\,j\right)\kappa$, and $\gamma=10^{-5}\kappa$.
Figure 9:  As in Fig. 8 (a) and (b) for the model “$Ccq$” with only qubits.

Finally, we analyze, in Figs. 8 and 9, the effect of random deviations of the
system parameters from the symmetric configurations identified in Secs. 2 and
3. The results show that, as expected, the steady state entanglement tends to
be reduced when the system is not symmetric. The reduction is more pronounced
for deviation in the values of the detunings than of the couplings, and
increases with the system size.

##  5 Conclusions

We have shown that it is possible to drive a quantum array into a pure steady
state featuring many entangled qubit pairs, by controlling the dissipative
dynamics of a single, central, element, either a cavity or a qubit. The steady
state is pure when dissipation acts only on the central element. When
additional decoherence is added on the other qubits, the stationary state is
an entangled mixed state, and its entanglement is large and robust provided
that the decay rate of the additional dephasing processes is much smaller than
that of the dissipative decay rate of the central element.

These models can be realized in a number of different physical situations with
atomic systems [42, 43, 44, 45, 46] and solid state nano-devices [47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59], allowing the realization of the chain
or of the star geometry.

An interesting related question is whether these approaches can be extended to
the preparation of more complex multipartite entangled qubit states, such as
the graph states which are the fundamental resource of measurement-based
quantum computers [60], in a way similar to what has been demonstrated for
bosonic modes in Ref. [22].

We acknowledge the support of PNRR MUR project PE0000023-NQSTI (Italy), and of
the European Union Horizon 2020 Programme for Research and Innovation through
the Project No. 862644 (FET Open QUARTET).

## References

  * [1] M. B. Plenio, S. F. Huelga, A. Beige, and P. L. Knight.  Cavity-loss-induced generation of entangled atoms  Phys. Rev. A, 59:2468, March 1999.  doi:10.1103/PhysRevA.59.2468.
  * [2] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller.  Preparation of entangled states by quantum Markov processes.  Phys. Rev. A, 78(4):042307, October 2008.  doi:10.1103/PhysRevA.78.042307.
  * [3] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller.  Quantum States and Phases in Driven Open Quantum Systems with Cold Atoms.  Nat. Phys., 4(11):878–883, November 2008.  doi:10.1038/nphys1073.
  * [4] F. Verstraete, M. M. Wolf, and J. I. Cirac.  Quantum computation and quantum-state engineering driven by dissipation.  Nat. Phys., 5(9):633–636, July 2009.  doi:10.1038/nphys1342.
  * [5] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler.  A Rydberg quantum simulator.  Nat. Phys., 6(5):382–388, May 2010.  doi:10.1038/nphys1614.
  * [6] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt.  An open-system quantum simulator with trapped ions.  Nature, 470(7335):486–491, February 2011.  doi:10.1038/nature09801.
  * [7] J. Cho, S. Bose, and M. S. Kim.  Optical Pumping into Many-Body Entanglement.  Phys. Rev. Lett., 106(2):020504, January 2011.  doi:10.1103/PhysRevLett.106.020504.
  * [8] K. Koga and N. Yamamoto.  Dissipation-induced pure Gaussian state.  Phys. Rev. A, 85(2):022103, February 2012.  doi:10.1103/PhysRevA.85.022103.
  * [9] G. Morigi, J. Eschner, C. Cormick, Y. Lin, D. Leibfried, and D. J. Wineland.  Dissipative Quantum Control of a Spin Chain.  Phys. Rev. Lett., 115(20):200502, November 2015.  doi:10.1103/PhysRevLett.115.200502.
  * [10] F. Reiter, D. Reeb, and A. S. Sørensen.  Scalable Dissipative Preparation of Many-Body Entanglement.  Phys. Rev. Lett., 117(4):040501, July 2016.  doi:10.1103/PhysRevLett.117.040501.
  * [11] S. Zippilli, M. Paternostro, G. Adesso, and F. Illuminati.  Entanglement Replication in Driven Dissipative Many-Body systems.  Phys. Rev. Lett., 110(4):040503, January 2013.  doi:10.1103/PhysRevLett.110.040503;  Erratum: Entanglement Replication in Driven Dissipative Many-Body Systems [Phys. Rev. Lett. 110, 040503 (2013)].  Phys. Rev. Lett., 111(16):169901, October 2013.  doi:10.1103/PhysRevLett.111.169901.
  * [12] S. Zippilli and F. Illuminati.  Non-Markovian dynamics and steady-state entanglement of cavity arrays in finite-bandwidth squeezed reservoirs.  Phys. Rev. A, 89(3):033803, March 2014.  doi:10.1103/PhysRevA.89.033803.
  * [13] S. Ma and M. J. Woolley.  Entangled pure steady states in harmonic chains with a two-mode squeezed reservoir.  J. Phys. A: Math. Theor., 52(32):325301, July 2019.  doi:10.1088/1751-8121/ab2ce7.
  * [14] P. Wendenbaum, T. Platini, and D. Karevski.  Entanglement replication via quantum repeated interactions.  Phys. Rev. A, 91(4):040303, April 2015.  doi:10.1103/PhysRevA.91.040303.
  * [15] A. Pocklington, Y.-X. Wang, Y. Yanay, and A. A. Clerk.  Stabilizing volume-law entangled states of fermions and qubits using local dissipation.  Phys. Rev. B, 105(14):L140301, April 2022.  doi:10.1103/PhysRevB.105.L140301.
  * [16] S. Zippilli, J. Li, and D. Vitali.  Steady-state nested entanglement structures in harmonic chains with single-site squeezing manipulation.  Phys. Rev. A, 92(3):032319, September 2015.  doi:10.1103/PhysRevA.92.032319.
  * [17] M. Asjad, S. Zippilli, and D. Vitali.  Mechanical Einstein-Podolsky-Rosen entanglement with a finite-bandwidth squeezed reservoir.  Phys. Rev. A, 93(6):062307, June 2016.  doi:10.1103/PhysRevA.93.062307.
  * [18] S. Ma, M. J. Woolley, I. R. Petersen, and N. Yamamoto.  Pure Gaussian states from quantum harmonic oscillator chains with a single local dissipative process.  J. Phys. A: Math. Theor., 50(13):135301, 2017.  doi:10.1088/1751-8121/aa5fbe.
  * [19] Y. Yanay and A. A. Clerk.  Reservoir engineering of bosonic lattices using chiral symmetry and localized dissipation.  Phys. Rev. A, 98(4):043615, October 2018.  doi:10.1103/PhysRevA.98.043615.
  * [20] Y. Yanay and A. A. Clerk.  Reservoir engineering with localized dissipation: Dynamics and prethermalization.  Phys. Rev. Research, 2(2):023177, May 2020.  doi:10.1103/PhysRevResearch.2.023177.
  * [21] Y. Yanay.  Algorithm for tailoring a quadratic lattice with a local squeezed reservoir to stabilize generic chiral states with nonlocal entanglement.  Phys. Rev. A, 102(3):032417, September 2020.  doi:10.1103/PhysRevA.102.032417.
  * [22] S. Zippilli and D. Vitali.  Dissipative Engineering of Gaussian Entangled States in Harmonic Lattices with a Single-Site Squeezed Reservoir.  Phys. Rev. Lett., 126(2):020402, January 2021.  doi:10.1103/PhysRevLett.126.020402.
  * [23] S. Ma, J. Zhang, X. Li, Y. Ren, J. Xie, M. Cao, and F. Li,  Coupling-modulation–mediated generation of stable entanglement of superconducting qubits via dissipation  EPL, 135:63001, 2021. doi:10.1209/0295-5075/ac2b5c.
  * [24] M. Asjad, S. Zippilli, and D. Vitali,  Suppression of Stokes scattering and improved optomechanical cooling with squeezed light  Phys. Rev. A, 94(5):051801, November 2016.  doi:10.1103/PhysRevA.94.051801.
  * [25] J. B. Clark, F. Lecocq, R. W. Simmonds, J. Aumentado, and J. D. Teufel.  Sideband cooling beyond the quantum backaction limit with squeezed light  Nature 541, 191–195 (2017).  doi:10.1038/nature20604.
  * [26] F. Acernese, et al. (Virgo Collaboration). Increasing the Astrophysical Reach of the Advanced Virgo Detector via the Application of Squeezed Vacuum States of Light. Phys. Rev. Lett., 123:231108, 2019. doi:10.1103/PhysRevLett.123.231108.
  * [27] M. Tse, H. Yu, et al.. Quantum-Enhanced Advanced LIGO Detectors in the Era of Gravitational-Wave Astronomy. Phys. Rev. Lett., 123:231107, 2019. doi:10.1103/PhysRevLett.123.231107.
  * [28] G. Vitagliano, A. Riera, and J. I. Latorre.  Volume-law scaling for the entanglement entropy in spin-1/2 chains.  New J. Phys., 12(11):113049, November 2010.  doi:10.1088/1367-2630/12/11/113049.
  * [29] G. Ramírez, J. Rodríguez-Laguna, and G. Sierra.  From conformal to volume law for the entanglement entropy in exponentially deformed critical spin 1/2 chains.  J. Stat. Mech., 2014(10):P10004, October 2014.  doi:10.1088/1742-5468/2014/10/P10004.
  * [30] C. M. Langlett, Z.-C. Yang, J. Wildeboer, A. V. Gorshkov, T. Iadecola, and S. Xu.  Rainbow scars: From area to volume law.  Phys. Rev. B, 105(6):L060301, February 2022.  doi:10.1103/PhysRevB.105.L060301.
  * [31] S. Zippilli, S. M. Giampaolo, and F. Illuminati.  Surface entanglement in quantum spin networks.  Phys. Rev. A, 87(4):042304, April 2013.  doi:10.1103/PhysRevA.87.042304.
  * [32] C. Di Franco, M. Paternostro, and M. S. Kim.  Nested entangled states for distributed quantum channels.  Phys. Rev. A, 77(2), February 2008.  doi:10.1103/PhysRevA.77.020303.
  * [33] B. Alkurtass, L. Banchi, and S. Bose.  Optimal quench for distance-independent entanglement and maximal block entropy.  Phys. Rev. A, 90(4):042304, October 2014.  doi:10.1103/PhysRevA.90.042304.
  * [34] I. Pitsios, L. Banchi, A. S. Rab, M. Bentivegna, D. Caprara, A. Crespi, N. Spagnolo, S. Bose, P. Mataloni, R. Osellame, and F. Sciarrino.  Photonic simulation of entanglement growth and engineering after a spin chain quench.  Nat Commun, 8(1):1569, November 2017.  doi:10.1038/s41467-017-01589-y.
  * [35] W. Cottrell, B. Freivogel, D. M. Hofman, and S. F. Lokhande.  How to build the thermofield double state.  J. High Energ. Phys., 2019(2):58, February 2019.  doi:10.1007/JHEP02(2019)058.
  * [36] A. R. Brown, H. Gharibyan, S. Leichenauer, H. W. Lin, S. Nezami, G. Salton, L. Susskind, B. Swingle, and M. Walter.  Quantum Gravity in the Lab: Teleportation by Size and Traversable Wormholes.  arXiv:1911.06314, February 2021.
  * [37] S. Dutta and N. R. Cooper.  Long-Range Coherence and Multiple Steady States in a Lossy Qubit Array.  Phys. Rev. Lett., 125(24):240404, December 2020.  doi:10.1103/PhysRevLett.125.240404.
  * [38] S. Dutta and N. R. Cooper.  Out-of-equilibrium steady states of a locally driven lossy qubit array.  Phys. Rev. Research, 3(1):L012016, February 2021.  doi:10.1103/PhysRevResearch.3.L012016.
  * [39] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki.  Quantum entanglement.  Rev. Mod. Phys., 81(2):865–942, June 2009.  doi:10.1103/RevModPhys.81.865.
  * [40] J. Onam González, L. A. Correa, G. Nocerino, J. P. Palao, D. Alonso, G. Adesso. Testing the Validity of the ’Local’ and ’Global’ GKLS Master Equations on an Exactly Solvable Model. Open Syst. Inf. Dyn., 24:1740010, 2017. doi:10.1142/S1230161217400108; P. P. Hofer, M. Perarnau-Llobet, L. David M. Miranda, G. Haack, R. Silva, J. Bohr Brask, N. Brunner. Markovian master equations for quantum thermal machines: local versus global approach. New J. Phys., 19:123037, 2017. doi:10.1088/1367-2630/aa964f; M. Cattaneo, G. L. Giorgi, S. Maniscalco, R. Zambrini. Local versus global master equation with common and separate baths: superiority of the global approach in partial secular approximation. New J. Phys., 21:113045, 2019. doi:10.1088/1367-2630/ab54ac; S. Scali, J. Anders, L. A. Correa. Local master equations bypass the secular approximation. Quantum, 5:451, 2021. doi:10.22331/q-2021-05-01-451.
  * [41] F. Benatti, R. Floreanini, and M. Piani.  Environment Induced Entanglement in Markovian Dissipative Dynamics.  Phys. Rev. Lett., 91(7):070402, 2003.  doi:10.1103/PhysRevLett.91.070402.; B. Kraus and J. I. Cirac.  Discrete Entanglement Distribution with Squeezed Light.  Phys. Rev. Lett., 92(1):013602, January 2004.  doi:10.1103/PhysRevLett.92.013602; M. Paternostro, W. Son, and M. S. Kim.  Complete Conditions for Entanglement Transfer.  Phys. Rev. Lett., 92(19):197901, May 2004.  doi:10.1103/PhysRevLett.92.197901. G. Adesso, S. Campbell, F. Illuminati, and M. Paternostro.  Controllable Gaussian-Qubit Interface for Extremal Quantum State Engineering.  Phys. Rev. Lett., 104(24):240501, June 2010.  doi:10.1103/PhysRevLett.104.240501.
  * [42] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe.  Observation of a Many-Body Dynamical Phase Transition with a 53-Qubit Quantum Simulator.  Nature, 551(7682):601–604, November 2017.  doi:10.1038/nature24654.
  * [43] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos.  Probing Rényi entanglement entropy via randomized measurements.  Science, 364(6437):260–263, April 2019.  doi:10.1126/science.aau4963.
  * [44] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller.  Self-verifying variational quantum simulation of lattice models.  Nature, 569(7756):355–360, May 2019.  doi:10.1038/s41586-019-1177-4.
  * [45] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage.  Trapped-Ion Quantum Computing: Progress and Challenges.  Applied Physics Reviews, 6(2):021314, June 2019.  doi:10.1063/1.5088164.
  * [46] M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti, T. Calarco, Z. Idziaszek, and P. S. Julienne.  Cold hybrid ion-atom systems.  Reviews of Modern Physics, 91:035001, July 2019.  doi:10.1103/RevModPhys.91.035001.
  * [47] S. Hacohen-Gourgy, V. V. Ramasesh, C. De Grandi, I. Siddiqi, and S. M. Girvin.  Cooling and Autonomous Feedback in a Bose-Hubbard Chain with Attractive Interactions.  Phys. Rev. Lett., 115(24):240501, December 2015.  doi:10.1103/PhysRevLett.115.240501.
  * [48] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck.  Observation of a Dissipative Phase Transition in a One-Dimensional Circuit QED Lattice.  Phys. Rev. X, 7(1):011016, February 2017.  doi:10.1103/PhysRevX.7.011016.
  * [49] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis.  Spectroscopic signatures of localization with interacting photons in superconducting qubits.  Science, 358(6367):1175–1179, December 2017.  doi:10.1126/science.aao1401.
  * [50] A. J. Kollár, M. Fitzpatrick, and A. A. Houck.  Hyperbolic lattices in circuit quantum electrodynamics.  Nature, 571(7763):45–50, July 2019.  doi:10.1038/s41586-019-1348-3.
  * [51] R. Ma, B. Saxberg, C. Owens, N. Leung, Y. Lu, J. Simon, and D. I. Schuster.  A dissipatively stabilized Mott insulator of photons.  Nature, 566(7742):51–57, February 2019.  doi:10.1038/s41586-019-0897-9.
  * [52] A. Smith, B. Jobst, A. G. Green, and F. Pollmann.  Crossing a topological phase transition with a quantum computer.  Phys. Rev. Research, 4(2):L022020, April 2022.  doi:10.1103/PhysRevResearch.4.L022020.
  * [53] E. Kim, X. Zhang, V. S. Ferreira, J. Banker, J. K. Iverson, A. Sipahigil, M. Bello, A. González-Tudela, M. Mirhosseini, and O. Painter.  Quantum Electrodynamics in a Topological Waveguide.  Phys. Rev. X, 11(1):011015, January 2021.  doi:10.1103/PhysRevX.11.011015.
  * [54] K. Le Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró.  Many-body quantum electrodynamics networks: Non-equilibrium condensed matter physics with light.  Comptes Rendus Physique, 17(8):808–835, October 2016.  doi:10.1016/j.crhy.2016.05.003.
  * [55] X. Gu, A. F. Kockum, A. Miranowicz, Y.-X. Liu, and F. Nori.  Microwave photonics with superconducting quantum circuits.  Physics Reports, 718–719:1–102, November 2017.  doi:10.1016/j.physrep.2017.10.002.
  * [56] G. Wendin.  Quantum information processing with superconducting circuits: A review.  Rep. Prog. Phys., 80(10):106001, 2017.  doi:10.1088/1361-6633/aa7e1a.
  * [57] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto.  Topological photonics.  Reviews of Modern Physics, 91:015006, January 2019.  doi:10.1103/RevModPhys.91.015006.
  * [58] I. Carusotto, A. A. Houck, A. J. Kollár, P. Roushan, D. I. Schuster, and J. Simon.  Photonic materials in circuit quantum electrodynamics.  Nat. Phys., 16(3):268–279, March 2020.  doi:10.1038/s41567-020-0815-y.
  * [59] S. A. Wilkinson and M. J. Hartmann.  Superconducting quantum many-body circuits for quantum simulation and computing.  Appl. Phys. Lett., 116(23):230501, June 2020.  doi:10.1063/5.0008202.
  * [60] M. Hein, W. Dür, J. Eisert, R. Raussendorf, M. Van den Nest, and H.-J. Briegel.  Entanglement in graph states and its applications.  Proc. Int. Sch. Phys. Enrico Fermi, 162(Quantum Computers, Algorithms and Chaos):115–218, 2006.  doi:10.3254/978-1-61499-018-5-115.