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

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.


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.

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 |ψ , and that of the cavity/ies |ϕ 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 − jth, and the steady state of the qubits can be expressed as 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), |± j is the eigenvalue of the Pauli operator σ (z) j , for the qubit j, with eigenvalue ±1, and n is the number of excitations of the squeezed reservoir.Moreover, χ j is a phase factor which depends on the specific model [as specified in Eq. ( 12)].We finally note that, by increasing n, the state of each pair in Eq. ( 1) tends to a Bell, maximally-entangled state, which is a central resource in many quantum information protocols [39].

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 "S cq" 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 "S q".

The master equation
In every case, the system dynamics is described by a master equation of the form with ξ ∈ {Ccq, S cq, Cq, S q}.The Lindblad operator L c describes dissipation via the squeezed reservoir of the central cavity with photonic annihilation operator b 0 (see Fig. 1), and it reads where 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].

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 ∝ b † j 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 ∝ b † j σ j +h.c., with σ j = σ (x) j − i σ (y) j /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 ∝ σ (x) j σ (x) k + σ (y) j σ (y) k = 2 σ † j σ k + h.c.. To be specific, the Hamiltonians corresponding to each of the four models are given by the expressions reported in Eqs. ( 5)-( 8) below, which describe the system in a reference frame rotating at the frequency of the central cavity ω 0 /2π.In particular -on the one hand -in the models "Ccq" and "S cq" 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 ∆ c, j from ω 0 .On the other hand, in the models "Cq" and "S q", composed of many qubits and a central cavity, the transition frequency of each qubit is detuned by ∆ q, j from ω 0 .Our analysis relies on the assumption that ω 0 is the dominant parameter in the studied systems.Consequently, the frequencies of the cavities and qubits (ω 0 + ∆ c, j and ω 0 + ∆ q, j , respectively) are orders of magnitude larger than the introduced coupling strengths, a common feature of quantumoptical 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. ( 1), 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)] 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 | j| 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: "S cq" [Fig. 1 (b)] 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)] 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: "S q" [Fig. 1 (d)] As in the "S cq" 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.

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 "S cq") approach a pure entangled state constituted of many twomode squeezed states.It can be expressed as where |0 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 "S cq"), 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 where tanh(r) = n/ ( n + 1).Instead, for the models "Cq" and "S q", U c is the single mode squeezing operator In particular the factor χ j , which appears in Eqs. ( 1) and (10), is given by for j 0. We also note that these operators [Eqs.(10) and (11)] realize the Bogoliubov transformation for all j, also for j = 0 with χ 0 = 1.Now, one can check that, in general, the product state between the state of the qubits given by Eq. ( 1), and that of the cavity/ies given by Eq. ( 9), with U c and χ j defined in Eqs. ( 10)-( 12), is a steady state for the four models.In other words, one finds that − i H ξ , |Ψ Ψ| + κ 2 L c |Ψ Ψ| = 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,ξ (this is non-zero only for the models "Ccq" and "S cq"), the one for the qubits alone H q,ξ (this is non-zero only for the models "Cq" and "S q") and the terms which describe the interactions between cavity/ies and qubits H c−q,ξ , such that we find that as demonstrated in Refs.[16,19,22], and because of the destructive quantum interference between transitions which involve the qubit states in the quantum superposition of Eq. (1).

Dynamics in the squeezed representation
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 ρ = U † c ρ U c , which is given by with the Lindblad term which describes dissipation of the zeroth mode being The transformed Hamiltonian in Eq. ( 19) can be written as where the Jaynes-Cummings interaction term H c−q,ξ = U † c H c−q,ξ U c may be expressed as Low-excitations eigenlevels (black horizontal lines) corresponding to the Hamiltonian H Ccq η c, j =0,g j =0 (in the squeezed representation) without interaction terms [see Eqs. ( 5), ( 21) and (22)] 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 "S cq" with N = 1).The purple wavy arrow illustrates the transfer of population due to the decay of the central cavity; 2 n + 1 is its orthogonal state.The parameters n j and n ≡ j n j indicate the Fock states |n −1 n 0 n 1 c of the cavities in the squeezed representation (i.e., the squeezed Fock states U c |n −1 n 0 n 1 c in the original representation).The values of λ tag the eigenvalues of H Ccq η 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 |ψ 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.
with τ j given by the collective qubit operator Now it is easy to verify that the transformed state (with |... = U † c |... ) fulfills the relation H c−q,ξ |Ψ = 0 [which is equivalent to Eq. ( 18)].In fact, on the one hand, b j |Ψ = 0 because the cavity(-ies) is (are) in the vacuum (in this representation) and, on the other hand, 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. ( 1) 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 "S cq" 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. ( 1) 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.

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 with other lowering and rising qubits operators σ c, j and σ † c, j (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).
Let us explicitly write down the equations for these models.The master equation has the form of Eq. ( 19) where also in this case we use the labels ξ ∈ {Ccq, S cq, Cq, S q} 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 → σ c, j .In other words k + σ (y) j σ (y) k , while the dotted black lines indicate interactions of the form ∝ σ (x) j σ (x) k − σ (y) j σ (y) k .In (c) and (d) the interactions between the central spin and the neighboring ones is of XY-type with anisotropic couplings ∝ g (x) σ (x) j σ (x) k + g (y) σ (y) j σ (y) k , while the black lines in (c) between the other spins indicate XXinteractions ∝ σ (x) j σ (x) k + σ (y) j σ (y) k . namely and where, as in Sec. 2, H q,ξ is non-zero only for ξ ∈ {Cq, S q} and H • c,ξ is non-zero only for ξ ∈ {Ccq, S cq}, with , , 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 "S cq", each qubit corresponding to a cavity of the previous model interacts with two qubits according to the Hamiltonians , for ξ ∈ {Ccq, S cq} .
Instead, in the models "Cq" and "S q", the Jaynes-Cummings terms result in the anistrotopic XY-interaction Hamiltonians where χ 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 where |− c indicates the state for all the qubits with lowering operator σ c, j , where each qubit is in the eigenstate of σ (z) c, j with eigenvalue −1.

Numerical results
The fact that the steady state Eq. ( 14) is in fact unique, for specific choices of the parameters ∆ c, j , η c, j , ∆ q, j , η 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 ‡, according to the equation ˙ ρ = L ρ, where the total Liouvillian superoperator is given by and where the additional noise on the qubits is described by the Lindblad term 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 ‡ The 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 κ. ∆ q,1 −0.193κ and g 1 = 0.36κ.In (b) and (e) ∆ q,1 −0.193κ, ∆ q,2 = ∆ q,1 + 0.05κ, g 1 = 0.36κ, and η q,1 = 0.362κ.In (c) and (f), the dimensions of the Hilbert spaces of all the cavities are truncated at the same Fock number n (max) 0 = n (max) ±1 , and ∆ c,1 −0.26κ, g 1 = 0.36κ.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 L as a function of ∆, see Fig. "S q", all the detunings and all the couplings are equal, such that ∆ q, j = ∆ q, j and g j = g j for all j, j ∈ {1, . . ., N}.In this highly symmetric case, any partition of the qubits in pairs may be used to construct a state as that of Eq. ( 1) 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 γ), 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.
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 (max) j (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 (max) j = 1 corresponding to models with only qubits of Sec. 3 and the largest values of n (max) j 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.
In Fig. 4 (a)-(c) we observe that in the ideal case (γ = 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. ( 1) 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 n.However, we notice that the dynamics involving cavities [corresponding to larger values of n (max) j ] is significantly faster than that of the models with qubits only [n (max) j = 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)].
Figs. 5 and 6 show that, on the one hand, at zero additional noise (γ = 0) maximum entanglement is achieved for larger values of n, as expected from Eq. ( 1), which approaches the product of many Bell, maximally entangled states.On the other hand, when the rate γ of the additional noise is finite, maximum entanglement is achieved at finite values of n.In fact, the larger n, the slower the dynamics, as illustrated by these figures [see the decay rates for different values of n in Fig. 5 (d)] and, as a consequence, if γ is too large noise and decoherence have enough time to spoil the generation of the steady state.
The fact that a larger 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 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 0.197κ.The coupling strength is g 1 = 0.36κ.The cavity Fock space is truncated at n (max) 0 = 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 (max) 0 = 1) and for larger values of n (max) 0 .(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 L identified by the vertical dashed line in plot (c).relations ∆ q, j = ∆ + 0.05κ( j − 1) and η q, j = g 1 + 0.002κ( j − 1).
In this figure, we also maximize the final concurrence as a function of ∆ only.Specifically, the value of ∆ is chosen in order to maximize the smallest real part of the eigenvalues of 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].
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.In (a) and (d) the solid (dashed) lines refer to the couple j = 2 (1).In (a)-(c) ∆ q,1 −0.13κ, ∆ q,2 = ∆ q,1 + 0.05κ, g 1 = 0.36κ, and η q,1 = 0.362κ.The Fock space of the cavity is truncated at n (max) 0 = 2.In (d)-(f) ∆ q,1 1.3κ, ∆ q,2 = ∆ q,1 + 0.05κ, g 1 = 0.36κ, and g 2 = 0.362κ.The Fock space of the cavity is truncated at n (max) 0 = 2.In (g)-(i) g 1 = 0.36κ and ∆ q,1 −0.17κ.The Fock spaces of the cavities are truncated at n (max) 0 = 2 and n (max) ±1 = 1.

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.
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].The Fock space of the central cavity is truncated to n (max) 0 = 2.The system parameters are g 1 = 0.36κ, η q, j = g 1 + 0.002κ( j − 1), and ∆ q, j = ∆ + 0.05κ( j − 1), where the value of ∆ 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.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 ∆ q, j we use ∆ (r)  q,± j = ±∆ q, j 1 + r (∆) ± j with r (∆) ± j random variables uniformly distributed in the range [−r max , r max ].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 η j we use g (r) ±1 = g 1 1 + r (g) ± j and η (r) ± j = η j 1 + r (η) ± j , respectively, with r (g) ±1 and r (η) ± j random variables uniformly distributed in the range [−r max , r max ].In each plot, for each value of r 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 = r (g) j = r (η) j = 0.In all plots ∆ j = (1.8+ 0.2 j) κ, g 1 = 0.3 κ, η j = (0.45 + 0.05 j) κ, and γ = 10 −5 κ.

Figure 1 .
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 "S cq".(c) Chain of qubits with a central cavity "Cq".(d) Star of qubits with a central cavity "S q".

Figure 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 ∝ σ (x) j σ (x)k + σ (y) j σ (y) k , while the dotted black lines indicate interactions of the form ∝ σ (x) j σ (x) k − σ (y) j σ (y) k .In (c) and (d) the interactions between the central spin and the neighboring ones is of XY-type with anisotropic couplings ∝ g (x) σ (x) j σ (x) k + g (y) σ (y) j σ (y) k , while the black lines in (c) between the other spins indicate XXinteractions ∝ σ (x) j σ (x) k + σ (y) j σ (y) k .

Figure 5 .
Figure 5. (a) Steady state concurrence for the model "Cq" with N = 1 (that is equal to "S q" with N = 1) as a function of the average number of excitations of the squeezed reservoir n, and for various values of the dephasing rate γ.(b) Time evolution of the concurrence for the highlighted (yellow) points in (a) corresponding to n = 0.5 (dark red), 1 (red), and 5 (salmon) with γ = 5 × 10 −4 κ.(c) Real part of the first two eigenvalues of L for the same points.The initial state is vacuum (in the squeezed representation) for the cavity and the eigenstate of σ (z) j with eigenvalue −1 for all the qubits.The vertical black, dashed line in (c) indicates the value of ∆ q,10.197κ.The coupling strength is g 1 = 0.36κ.The cavity Fock space is truncated at n(max)

Figure 9 .
Figure 9.As in Fig. 8 (a) and (b) for the model "Ccq" with only qubits.