Dissipative Bose-Einstein condensation in contact with a thermal reservoir

We investigate the real-time dynamics of open quantum spin-$1/2$ or hardcore boson systems on a spatial lattice, which are governed by a Markovian quantum master equation. We derive general conditions under which the hierarchy of correlation functions closes such that their time evolution can be computed semi-analytically. Expanding our previous work [Phys. Rev. A 93, 021602 (2016)] we demonstrate the universality of a purely dissipative quantum Markov process that drives the system of spin-$1/2$ particles into a totally symmetric superposition state, corresponding to a Bose-Einstein condensate of hardcore bosons. In particular, we show that the finite-size scaling behavior of the dissipative gap is independent of the chosen boundary conditions and the underlying lattice structure. In addition, we consider the effect of a uniform magnetic field as well as a coupling to a thermal bath to investigate the susceptibility of the engineered dissipative process to unitary and nonunitary perturbations. We establish the nonequilibrium steady-state phase diagram as a function of temperature and dissipative coupling strength. For a small number of particles $N$, we identify a parameter region in which the engineered symmetrizing dissipative process performs robustly, while in the thermodynamic limit $N\rightarrow \infty$, the coupling to the thermal bath destroys any long-range order.


I. INTRODUCTION
With the advent of ultracold gases and trapped ions as tunable quantum simulators, experiments are now in the position to investigate the real-time evolution of quantum manybody systems directly with engineered model Hamiltonians [1][2][3]. Recent years have seen tremendous progress that promises new insights at the intersection of condensed matter physics, high energy physics, and beyond [4,5]. While experiments have been very successful in probing the dynamics of quantum many-body systems, it is fair to say that to date our theoretical understanding remains incomplete. Thus, many important problems regarding nonequilibrium quantum dynamics as, e.g., the characterization of steady states or the calculation of transport properties far from equilibrium are still beyond our reach.
There are several reasons for this disparity. While exact diagonalization provides us with rigorous results it is limited to systems consisting of a small number of particles. The density matrix renormalization group (DMRG) [6,7] allows for the efficient simulation of real-time dynamics for a large variety of one-dimensional gapped quantum many-body systems [8][9][10][11]. Nevertheless, its applicability remains largely limited to short times due to the growth of entanglement. Quantum Monte Carlo (QMC) methods are not limited by macroscopic quantum correlations, however, their application to nonequilibrium dynamics is hindered by a severe complex phase problem.
In light of these limitations, it is quite remarkable that a class of open quantum manybody models [12] can be solved without approximations. This is of particular interest since open quantum system with engineered couplings to an environment [13] have been proposed recently for the preparation of quantum states [14][15][16][17][18], quantum simulation [19][20][21], as well as quantum computing [14,[22][23][24]. Typically, the nonequilibrium steady state (NESS) of the engineered dissipative process, which defines a unique fixed point for the dynamics, is known by construction; 1 examples include condensates of bosons or η states of fermions [25,26], condensates of hardcore bosons or quantum spins [27], d-wave pairing states of fermions [28,29], and various topologically ordered states [30][31][32][33]. So far, for most of these systems the real-time evolution leading to the final state is less well understood. It is known, however, that the dissipative contributions to the dynamics might in fact close the hierarchy of n-particle correlation functions such that a semi-analytic solution of interacting quantum systems becomes feasible. This observation has been exploited in a number of works [34][35][36][37][38][39][40] and general closure conditions for the hierarchy of correlation functions have been derived for bosonic and fermionic models [41].
In recent work [42], we investigated a protocol for dissipative generation of a Bose-Einstein condensate (BEC) [25][26][27]. In particular, we considered a macroscopic system of N spin-1/2 particles, where the time evolution of the reduced density matrix ρ is governed by a Markovian quantum master equation [12,43,44] λ = 1, . . . , d 2 N − 1, where d N = 2 −N denotes the dimension of the N-particle Hilbert space. The Lindbladian L = −i[H, · ] + λ L λ is a linear map on the set of mixed states and the commutator part defines the system Hamiltonian H, which is Hermitian. The dissipative coupling to the environment is described by the non-commutator part λ L λ of L. Written in Lindblad form, the action of its constituents L λ on ρ is expressed in terms of jump where the rate parameters γ λ ≥ 0 characterize the relative strength of the dissipative couplings. Eq. (1) represents an effective description of the system where the environment has been integrated out, yielding operators √ γ λ L λ that are expressed in terms of the local spin degrees of freedom s α x , where α = 1, 2, 3 denotes the spin index and x labels distinct particles. Specifically, in Ref. [42] only a single class of non-Hermitian jump operators was considered, which acts on pairs of particles, i.e., L xy = 1 In the absence of a Hamiltonian contribution (H = 0), we were able to show that this process yields a closed hierarchy of correlations, thereby allowing us to solve for the full real-time evolution of spin-spin correlation functions. Thus, the dynamics of dissipative Bose-Einstein condensation could be followed explicitly by studying the relaxation towards the final state ρ NESS = |BEC BEC|, satisfying Lρ NESS = 0.
Here, we go beyond this first exploratory work and establish the generality of our approach, which allows us to investigate the susceptibility of the considered quantum Markov process to unitary and nonunitary perturbations. Specifically, we consider operators L x that act locally on particles at site x, as well as bilocal jump operators L xy that act on pairs of particles (x, y), both in the presence and absence of a uniform magnetic field. While multilocal operators L x 1 x 2 ···xn , with n > 2, are conceivable and have been studied theoretically, e.g., in the context of steady states with nontrivial topology [45][46][47][48], local and bilocal operators appear to be sufficient to describe the phenomenology of engineered dissipative quantum spin-1/2 systems for BEC generation [27,49].
The outline of this paper is as follows: In Sec. II we derive the closure conditions for hierarchies of correlation functions of quantum spin systems in the s = 1/2 representation in the case of open Markovian dynamics governed by local and bilocal jump operators. In Sec. III we discuss the dynamics of an engineered dissipative process that drives the system into a mixture of totally symmetric superposition states. We describe the growth of longrange correlations in real time and consider the finite-size scaling of the dissipative gap. In Sec. IV we investigate the same dissipative process in the presence of thermal noise and a uniform magnetic field. Our results are summarized in the nonequilibrium steady-state phase diagram for the coupled model. We conclude with an outlook on applications and possible further developments of this work.

II. CLOSED HIERARCHIES FOR CORRELATION FUNCTIONS
Throughout this work we consider quantum spin-1/2 systems on a regular d-dimensional lattice. The spin degrees of freedom s α x = 1 2 σ α x are expressed in terms of Pauli matrices defined locally at site x, while the total spin operator is given by S α = x s α x . The assumption of an underlying regular lattice structure is not essential, however, it allows us to exploit the symmetries of the lattice to simplify the problem and to solve for the dynamics of correlation functions of macroscopic quantum many-body systems.
In the following, we consider the dynamics of n-point correlation functions O z 1 z 2 ···zn = tr {ρO z 1 z 2 ···zn }, in particular, products of spin operators O z 1 z 2 ···zn = 1≤i≤n s z i (where spin indices α i are omitted). Using Eq. (1) we derive the equation of motion where we have introduced the dual (or adjoint) map L * = i[H, · ] + λ L * λ corresponding to the Lindbladian L. In general, Eq. (3) does not close, i.e., the commutators on the right hand side typically induce operators of higher order. This leads to an infinite hierarchy of correlation functions which cannot be solved without truncating the coupled set of equations.
Here, we seek conditions under which the opposite is true. That is, we derive conditions under which the hierarchy of correlation functions closes for a purely dissipative quantum Markov process (H = 0) defined in terms of local (Sec. II A) or bilocal jump operators (Sec. II B). We discuss possible solutions to these conditions and return to the effect of unitary perturbations in Sec. II C.

A. Local jump operators
We start with the discussion of local operators √ γ x L x (γ x > 0) to illustrate the derivation of closure conditions. Any local jump operator can be expressed in the form where the coefficients l 0 and l α , α = 1, 2, 3, are defined up to an overall phase factor, which we may use to set l 0 ∈ Ê and l α ∈ . Since L x acts locally and spin operators at different sites commute, we observe that the commutator terms in Eq. (3) only contribute when x = z i , i = 1, . . . , n (cf. Fig. 1a). Thus, the equation of motion for n-point correlation functions can be written as with the dual superoperator L * which maps the spin s α z at site z = x onto the complete basis of the local spin-1/2 operator algebra, {½, s 1 z , s 2 z , s 3 z }, and to zero otherwise (z = x). We may express its action in the following form: where the additional factor 1/γ x is introduced to render the right hand side dimensionless.
The coefficients m α and m αβ are obtained by using the parametrization of the local jump operator [cf. Eq. (4)] to evaluate the action of L * x ; we find and the bar, e.g.,l α , denotes complex conjugation. We observe that the equation of motion closes, i.e., the equations of motion for n-point correlation functions do not couple to higher order correlations, independent of the particular form of the local jump operators L x . This does not hold true for generic jump operators that act on multiple sites, as we illustrate in the next section.

B. Bilocal jump operators
A large number of suggested protocols for dissipative state generation rely on Lindblad dynamics induced by jump operators that act on pairs of particles (see, e.g. Ref. [32]). It is therefore interesting to ask under what conditions the real-time dynamics can be solved either exactly or semi-analytically. Here, we restrict ourselves to bilocal operators √ γ xy L xy (γ xy > 0) that act isotropically on adjacent lattice sites x and y (which we denote by x, y ).
Although the assumptions of isotropy and nearest-neighbor couplings are both reasonable, they are not necessary for the following derivations. They serve merely to make the following arguments more transparent and to simplify the discussion.
For isotropic systems, any jump operator L xy is either symmetric (L xy = L yx ) or antisymmetric (L xy = −L yx ) under the interchange of x and y. By making use of this fact, we may reduce the number of coefficients that parametrize any bilocal operator L xy and thereby provide explicit closure conditions for such operators.

Symmetric jump operators
We employ the following parametrization for symmetric operators where we have defined the symmetrization of indices by l (αβ) = 1 2 (l αβ + l βα ) and the coefficients can be chosen as l 0 ∈ R and l α , l αβ ∈ C. The factor of i is introduced for later convenience. Specific examples for symmetric jump operators, whose purely dissipative dynamics has been investigated recently via a novel Monte Carlo method [50][51][52], are the singlet projection operator and triplet projection operator respectively.

Antisymmetric jump operators
In the case of antisymmetric operators, we may choose the following parametrization where l [αβ] = 1 2 (l αβ − l βα ) and the parameters l α and l αβ are complex valued in general but due to the phase ambiguity, we may choose any one of these parameters to be real. The non-Hermitian operator which provides a mechanism for dissipative cooling into a BEC [25][26][27]42] falls into this class of operators.

Closure conditions
Before we go on to discuss the closure conditions for n-point spin correlation functions, we derive the equation of motion for the local magnetization O z = tr{ρs z } and investigate under which conditions it decouples from higher order correlation functions.
Evaluating Eq. (3) for n = 1 we find that the commutator terms on the right hand side contribute only when jump operators L yz are attached to the site z (cf. Fig. 1b). In fact, these contributions have to be summed over all nearest-neighbor pairs y, z , so that The dual map L * xy is defined in terms of bilocal operators, i.e., L * xy = γ xy /2 L † xy · , L xy + L † xy , · L xy and its action on single spin operators can be expressed in the form Note that L * xy is invariant under the interchange of sites x and y. Clearly, the equation of motion for the one-point function closes only when the contributions ∼ s x s y in Eq. (14) vanish. Hence, to derive the closure conditions for bilocal operators L xy , we require that m αβγ = 0, which yields a set of 27 equations. Due to the distinct parametrization of the symmetric and antisymmetric jump operators, we discuss both cases separately: 1) We obtain the following set of closure conditions for symmetric jump operators [cf. Eq.
where we have introduced the coefficients a α = l αα , with α = 1, 2, 3, as well as b 1 = l (23) , (31) , and b 3 = l (12) . From these relations we see that a 10-parameter family of symmetric jump operators can be found (with l 0 , l α , a α , b α ∈ Ê), for which the above conditions are fulfilled and the local magnetization decouples from the higher order npoint correlation functions (n ≥ 2). Both the singlet and triplet projection operators P s xy and P tr xy belong to this class of operators. 2) For the antisymmetric jump operators [cf. Eq. (11)] we obtain where k 1 = l [23] , k 2 = l [31] , and k 3 = l [12] . We find a 6-parameter family of antisymmetric operators L xy (with l α , k α ∈ R) for which the equation of motion of the one-point function closes. Clearly, the non-Hermitian operator Q xy belongs to this class.
It remains to be shown that Eqs. (15a) -(15g) and Eq. (16) are sufficient to close the hierarchy for arbitrary n-point correlation functions. In the following we demonstrate this explicitly for two-point functions, but similar arguments also apply to correlation functions of higher order. The equation of motion for O z 1 z 2 = tr{ρs z 1 s z 2 } is given by where δ z 1 ,z 2 is equal to one only if z 1 and z 2 are nearest-neighbor sites. Note that the contributions on the right hand side can be simplified by using, e.g., Imposing the closure conditions for the local magnetization, i.e., Eqs. Hamiltonian that takes the following form with real-valued coefficients h α , J αβ = J βα . We evaluate the first term on the right hand side in Eq. (3), which yields Calculating the ensemble average, we see that the second term on the right hand side will induce a dependence on higher order correlation functions. Thus, to close the hierarchy we require that J αβ = 0. Accordingly, we will only allow for local Hamiltonian contributions in the following.

III. DISSIPATIVE COOLING INTO A BEC
In previous work [42], we studied the purely dissipative dynamics of a many-body quantum spin system on a hypercubic d-dimensional lattice with periodic boundary conditions.
The considered dissipative process is distinguished by non-Hermitian jump operators that lead to the growth of correlations and the generation of macroscopic order in the final state.
We discussed the dependence of the dissipative gap on the system size and found a novel nontrivial scaling behavior as a function of the system size N = L d , where L denotes the linear extent. In this section, we briefly summarize our main results and comment on the universality of this process with respect to the chosen lattice structure and boundary conditions of the problem. In Sec. IV, we augment this investigation by including the effect of thermal noise and studying the stability of the dissipative process.
A. Non-Hermitian jump operator and nonequilibrium steady state We consider a quantum Markov system of spin-1/2 particles that is driven uniformly on the lattice by the application of non-Hermitian operators √ γQ xy , where which acts on adjacent lattice sites x, y . The associated dissipative coupling, denoted by γ, is independent of the particular pair of particles x, y . The operator Q xy maps any two-particle spin-singlet state to the spin triplet, while conserving the total spin projection This state is unique for a given initial state and corresponds to an ensemble of totally symmetric superposition states [25]. By virtue of the quantum spin-1/2 to hardcore boson mapping [53], this dissipative process can also be seen as a symmetric delocalization of hardcore bosons over adjacent sites, with a BEC of hardcore bosons as the resulting final state.

B. Evolution equations for correlation functions
In Sec. II B we have seen that the hierarchy of correlation functions closes for the non-Hermitian jump operator Q xy . Here, we derive the explicit form for the equations of motion of one-and two-point functions when the dynamics is governed only by the dissipative process √ γQ xy on nearest neighbor sites, neglecting the effect of Hamiltonian contributions.

Local magnetization
We begin with the time evolution equation for the local magnetization S α x = tr {ρs α x }. To derive the equation of motion, we use the relation and sum over all nearest-neighbor pairs x, y , keeping x fixed. Calculating the trace over the density matrix ρ (using Eq. (3)) we obtain the diffusion equation The dimensionless time variable τ = γ t is introduced for later convenience, as well as the discretized Laplacian ∆ x , i.e., ∆ x δ xy = d µ=1 δ (x+μ)y − 2δ xy + δ (x−μ)y , where the lattice constant is set to one, andμ denotes the unit vector in the µ direction on the regular spatial lattice.

Spin-spin correlation functions
To follow the growth of long-range correlations we consider the ensemble average of bilocal spin operators s α x s β y . In particular, we consider C xy = s + x s − y +s − x s + y and D xy = s 3 x s 3 y , for which we define the expectation values C xy = tr {ρC xy } and D xy = tr {ρD xy }. Note that similar to the local magnetization, we will use calligraphic fonts to denote ensemble-averaged quantities in the following. The adjoint map L * xz corresponding to the non-Hermitian jump operator Q xz maps the operators C xy and D xy (x = y) to while the diagonal elements C xx = 4D xx = ½ are mapped to zero. The time evolution equation for the expectation values C xy and D xy are obtained by summing over all nearestneighbor pairs of x and y, and averaging over ρ [cf. Eq. (17)]. For x = y, we obtain d dτ while the diagonal terms C xx = 4D xx = 1 are constant in time. Given initial data for the two-point correlation functions, the above first-order system of differential equations can be solved explicitly via matrix diagonalization, where M denotes the linear differential operator that is determined by Eq. (26a) and (26b). Its solutions can be expressed in terms of a linear combination of exponential functions, whose characteristic decay rates are determined by the nonvanishing eigenvalues of M (see Sec. III C).
We exploit spatial translation invariance in the following to explicitly solve the linear system of equations (27). Accordingly, we may characterize the time evolution of spin-spin correlation functions in momentum space, e.g., with p µ = 2πn µ /L and n µ ∈ {0, . . . , L − 1}. The zero-momentum component C p=0 corresponds to the condensate fraction, which indicates the buildup of long-range correlations in the system.

C. Dynamics of purely dissipative cooling
As an example, we illustrate the time evolution starting from an incoherent thermal ensemble at infinite temperature also Fig. 2 in [42].
For this initial state, the spin-spin correlation functions are easily evaluated, i.e., indicating the absence of long-range correlations in the system. Also the final state can be found explicitly, which is expressed in terms of totally symmetric Dicke states, |D(N, n) = |N/2, −N/2 + n , which are simultaneous eigenstates of α S α 2 and S 3 , respectively. Using Eq. (31) we determine the asymptotic values for the correlation functions To solve for the full real-time evolution we diagonalize Eq. (27) numerically. Our results are shown in Fig. 2a, where we display the time evolution of the condensate fraction C p=0 (τ ).
The timescale on which the asymptotic regime is reached, i.e., C p=0 > C ∞ p=0 e −ǫ , 0 < ǫ ≪ 1, depends strongly on dimension. We may understand this behavior by inspecting the spectrum of the linear operator M [see Eq. (27)]. The corresponding eigenvalues λ ∈ Spec M take only nonpositive values and the mode corresponding to the eigenvalue which is denoted as dissipative gap in the following, dominates the asymptotic approach towards the final steady state. Interestingly, we find by numerically diagonalizing the full Lindbladian L on small system sizes that the value of the smallest negative eigenvalue of the linear operator M is identical to the dissipative gap of the Lindbladian L. It seems that the spin-spin correlation functions capture the slowest mode of the full Lindbladian L for the system under consideration.
For the purely dissipative process governed by the jump operators Q xy , the coupling γ can be scaled out, so that the total number of particles N = L d provides the only scale in the problem. Thus, the dissipative gap is a function of the system size and we may inquire about its asymptotic scaling properties, i.e., ∆ ∼ L −z , for sufficiently large L. We remark that the finite-size scaling of ∆ has been studied in detail for various one-dimensional bosonic and fermionic systems [34, 35, 37-39, 41, 54-56, 59]. However, here we observe that the considered purely dissipative quantum spin system exhibits a highly nontrivial scaling behavior that depends strongly on dimension [42]. The asymptotic finite-size scaling is given which is illustrated in Fig. 2b for the d = 1, d = 2 (square), and d = 3 (primitive cubic) lattice geometries.

D. Universality of the dissipative cooling process
The observed nontrivial finite size scaling of the dissipative gap for different dimensions d is remarkable. This concerns in particular the logarithmic correction in d = 2. One might speculate whether this behavior is due to the presence of topological defects (e.g., vortices) in the system. The deeper mathematical reason for this scaling behavior and its relation to the properties of M has not been fully elucidated yet. Consequently, we checked that this scaling behavior is not an artifact of the lattice structure or the boundary conditions of the problem.
In Fig. 3a, we display the scaling of ∆ for d = 2 dimensions upon changing the lattice geometry for fixed (periodic) boundary conditions. Specifically, we consider three different regular tilings of the plane, corresponding to a square ( ), triangular (▽), and honeycomb ( ) lattice geometry, respectively, and clearly observe the ∆ −1 ∼ N log N scaling for large N in all cases. However, note that the numerical values of ∆(N) differ. Empirically, we find that ∆ ▽ > ∆ > ∆ (see Fig. 3a), which can be attributed to the decrease of coordination number, i.e., n c,▽ > n c, > n c, .
Finally, we also investigated how the choice of boundary condition influences the finitesize scaling of the dissipative gap ∆. In Fig. 3b, we consider a d = 2 system on a square lattice with either open or periodic boundary conditions. Note that we are restricted to comparatively small lattice sizes when translation invariance is not imposed on the system. Nevertheless, we clearly observe that the above scaling is not altered by the choice of boundary conditions and we find ∆ pbc > ∆ obc . Summarizing, our results indicate that the dynamics of the purely dissipative cooling process for d = 2 is rather insensitive to the lattice structure or boundary conditions, and we expect a similar behavior for d = 1 and d = 3.

IV. DISSIPATIVE DYNAMICS IN THE PRESENCE OF COMPETING UNI-TARY AND NONUNITARY PROCESSES
Here, we supplement the purely dissipative cooling process considered in Sec. III by competing unitary and nonunitary processes. We restrict ourselves to Hamiltonians that include only a coupling to a uniform magnetic field, which is a necessary requirement to close the hierarchy of correlation functions. Moreover, we study the effect of thermally induced spin flips via an additional nonunitary process.

A. Competing unitary dynamics in the presence of a magnetic field
We consider the effect of an external field that points in the 1-direction, as described by where the operators E xy = s 2 x s 3 y + s 3 x s 2 y and F xy = s 2 x s 2 y have been introduced, in addition to C xy and D xy to close the set of evolution equations. The dissipative contributions that arise from the action of the adjoint map L * xz are proportional to Proceeding along the same lines as in the previous section, we obtain the full set of evolution where η = h/γ is the dimensionless magnetic field variable, rescaled by the dissipative coupling γ. The diagonal elements are constant in time, C xx = 4D xx = 4F xx = 1, E xx = 0, and the system of first-order differential equations (38a) -(38d) can be solved via diagonalization of the corresponding linear differential operator M η . We assume an incoherent thermal ensemble at infinite temperature as the initial state, such that C xy (0) = 4D xy (0) = 4F xy (0) = δ xy and E xy (0) = 0.
We first investigate the asymptotic behavior of the system by studying the spectrum of the linear differential operator M η , which is defined by the system of Eqs. (38a) -(38d). are rather seen to exhibit oscillations with frequency 2η around an average asymptotic value -the system is asymptotically characterized by a limit cycle. We determine the late-time behavior analytically for τ ≫ 1/∆ from Eqs. (38a) -(38d): The gray region indicates the time region for which τ < 1/∆.
Here, ϕ is an irrelevant phase offset and the function g(η, N) describes the oscillation amplitude at late times. In Fig. 4 we display the time evolution of the various two-point correlation functions, which clearly exhibits the oscillatory behavior at late times. Moreover, we numerically determine g(η, N) for which we observe a monotonic decay with increasing magnetic field strength η, where g(η → 0, N) = 1/16 and g(η → ∞, N) = (8Nη) −1 .
Finally, we display the time evolution of the condensate fraction C p=0 (τ ) for different values of η in Fig. 5. According to (39a), its late-time behavior is determined by In comparison to the purely dissipative cooling dynamics (η = 0, cf. Sec. III) for which the steady-state condensate fraction C ∞ p=0 = (N + 1)/(2N) is approached, we find that the presence of a nonvanishing magnetic field substantially decreases the late-time average value

Model two-spin system
To understand the main features of the dissipative real-time dynamics of Eqs. (39a) and (40), we may consider a model system consisting of two particles for which the time evolution of the density matrix can be easily calculated explicitly. To this end, we provide the initial density matrix with the triplet states |t + = |↑↑ , |t − = |↓↓ , |t 0 = (|↑↓ + |↓↑ )/ √ 2, as well as the singlet state |s = (|↑↓ − |↓↑ )/ √ 2. For vanishing magnetic field, the time-dependent density matrix is given by i.e., the singlet state |s s| is mapped to the triplet state |t 0 t 0 | while the other components remain invariant under the time evolution. This picture changes in the presence of a nonvanishing magnetic field. While the dissipative process of Sec. III still eliminates the singlet component from the ensemble, the magnetic field results in a mixing of the three triplet states where the time-dependent coefficients are given by Defining the two-particle operator C = s + ⊗ s − + s − ⊗ s + , we evaluate C |t + = C |t − = 0, C |t 0 = |t 0 , and C |s = − |s . Thus, it is sufficient to consider only the diagonal elements of the density matrix, denoted by ρ diag , to calculate the expectation value C = tr{ρC}. In fact, at late times, the time-averaged diagonal entries ρ diag are given by This result should be compared to the nonequilibrium steady-state density matrix in the absence of an external field Thus, the magnetic field generates an average spin rotation, which decreases the |t 0 t 0 | contribution when compared to the final state of the driven, purely dissipative system. As a consequence, the late-time average value lim τ →∞ 1 τ τ ∆ −1 dτ ′ C η (τ ′ ) = 3/8 is smaller than the asymptotic value for the purely dissipative cooling dynamics C η=0 = 1/2. It therefore appears that the phenomenology of the macroscopic N-particle system in the presence of a magnetic field in the 1-direction is essentially captured by the corresponding two-spin model system.

B. Competing thermal noise
Here, we introduce the effect of an external magnetic field that points in the 3-direction where we assume that h > 0. In contrast to H 1 , the Hamiltonian H 3 commutes with Q xy and therefore it does not lead to an additional coupling to two-point operators (as we encountered in the previous section). Here, however, we allow for local spin flips via additional jump processes that are accounted for by local operators L ± x = s ± x [56][57][58][59]. In general, we may assign independent interaction rates γ ± x to both processes L ± x . Assuming a thermal occupation of the bath, however, the spin flip rates are related via the Boltzmann where T denotes the bath temperature and we assume spatial homogeneity (γ ± x = γ ± ). This relation does not set the overall interaction rate which we denote by κ. We assign the following values to the ratios where n T ≡ e 2h/T − 1 −1 is the thermal occupation number. The equations of motion for correlation functions (3) receive additional contributions from the spin flip processes L ± x : while C xx = 4D xx = 1. The behavior of this system, which is driven by the operators √ γ ± L ± x and √ γQ xy can be fully characterized in terms of two independent parameters -the ratio of couplings κ/γ and the effective temperature T /h. To understand the relevant modes that determine the late-time behavior, we inspect the spectrum of the linear differential operator M T as defined by the system of linear equations (52a) -(52c). Again, we consider an incoherent infinitetemperature ensemble as initial state, for which C xy (0) = 4D xy (0) = δ xy and S 3 x (0) = 0. The spectrum of M T is real and nonpositive so that the asymptotic behavior is governed by the dissipative gap ∆ T . It is independent of the system size which is in stark contrast to purely dissipative cooling dynamics (κ = 0), for which the dissipative gap ∆ exhibits a nontrivial finite-size scaling. This means that the presence of a thermal bath dominates the asymptotic behavior of the system. More specifically, for the purely dissipative process governed by the jump operators Q xy , the three largest eigenvalues are given by λ 1 = λ 2 = 0 and λ 3 = −∆. The presence of thermal noise, however, lifts the degeneracy of the zero eigenvalues such that λ 1 = 0, λ 2 = −∆ T , and λ 3 < λ 2 . In Fig. 6 we compare the spectrum of a dissipative system in both cases -with and without thermal noise -and we display the scaling of the largest eigenvalues as a function of N. We clearly observe that ∆ T is independent of the system size when compared to ∆ (with κ = 0), which exhibits the nontrivial ∆ −1 ∼ N log N scaling in d = 2 dimensions.
Owing to the fact that L ± x does not conserve spin, [s 3 x , L ± x ] = 0, we obtain a nonvanishing value both for the asymptotic magnetization S 3 x and for the two-point function D xy : While these values can be calculated easily by hand, to determine C ∞ xy we need to invert the linear differential operator corresponding to the linear system Eqs. (52a) -(52c), which for large system sizes, can only be done numerically. We observe that the correlation function C ∞ xy (κ = 0) decays with the separation |x − y| in the presence of a thermal coupling, as shown in Fig. 7. This is in contrast to the dissipative process governed by the operator Q xy , for which the asymptotic value is given by C ∞ xy = 1/2, for x = y. We conclude that the thermally induced spin flips counteract the cooling process and destroy the long-range order in the system, thereby introducing a correlation length ξ ≪ L. The asymptotic value of the zero mode C ∞ p=0 strongly depends on the parameters γ/κ and T /h, as well as the particle number N. Setting the number of particles to N = 4096, we calculate the nonequilibrium steady-state phase diagram, which is shown in Fig. 8a. In this case, we find that the thermal noise destroys the long-range order more or less completely in the range γ/κ O(10 3 ), independent of the temperature T . This means that the coupling of the cooling process γ needs to be larger than that of the thermal bath κ by several orders of magnitude in order to generate a macroscopically ordered state. On the other hand, for large values of γ/κ, we observe an intriguing dependence of C ∞ p=0 on the temperature T [cf. Fig. 8a]. That is, for fixed γ/κ and finite N, we find a non-monotonous behavior of C ∞ p=0 (T ), where in the limiting cases lim T →0 C ∞ p=0 (T ) = lim T →∞ C ∞ p=0 (T ) = 1/N. In fact, these values corresponds to the lower bound for the zero mode C p=0 . This can be understood in the following way: In the limit T → 0, spin-flip operations s + x are forbidden, such that the action of s − x on all possible sites aligns all spins along the negative 3-direction, cf. Eq. (54), hence destroying any long-range order in the (1, 2)-plane. In the opposite limit T → ∞, we have γ + = γ − → ∞ such that the spin-flip operators s + x and s − x completely dominate the dynamics and, hence, prohibit any long-range order. Between these two limits, for any value of γ/κ, we observe that there exists an optimal value T opt for which C ∞ p=0 takes its maximum value [cf. Fig. 8b], also indicated by the continuous line in Fig. 8a. Thus, we may find a distinct temperature region T 1 < T < T 2 in which the engineered dissipation performs robustly even in the presence of thermal noise.
This region depends on the ratio of couplings γ/κ and the system size N. As the number of particles N is increased, the nearest-neighbor symmetrizing action Q xy finds it more and more difficult to compete against the thermal coupling that acts locally. Thus, we expect that in the thermodynamic limit N → ∞, single spin flips eventually destroy any long-range order and therefore dominate for any finite value of γ/κ. In contrast, for a small number of particles N the value of γ/κ that is required to generate long-range order decreases as well [cf. Fig. 8b]. In view of experimental realizations [27], our findings indicate that thermal fluctuations are not too prohibitive for generating long-range order via engineered dissipation, at least for not too large systems.

V. CONCLUSIONS
In this work we have studied the time evolution of correlations in the context of an open Markovian quantum many-body system, which was originally proposed for the dissipative cooling into a BEC [25][26][27]. In a previous publication [42] we showed that the corresponding purely dissipative process governed by the single non-Hermitian quantum jump operator Q xy = 1 2 s + x + s + y s − x − s − y allows for a semi-analytic solution of spin-spin correlation functions. Here, we have extended these results by studying the universality of the dissipative process.
We have established that the novel finite-size scaling behavior of the dissipative gap is in fact insensitive to the choice of lattice discretization as, e.g., provided by the lattice geometry or boundary conditions. Furthermore, we have studied the stability of the dissipative cooling process with respect to unitary and nonunitary perturbations. To this end, we derived conditions under which additional perturbations can be considered within the framework of a closed hierarchy of correlation functions, thereby admitting a closed analytic solution to the nonequilibrium dynamics. In particular, we allowed for the presence of a uniform magnetic field, as well as a coupling to a thermal bath κ, while the system is driven by the dissipative cooling process with a uniform rate γ. We calculated the nonequilibrium steady-state phase diagram and found that for any finite particle number N, above a certain threshold value for the coupling ratio γ/κ, the system allows for a final state ρ NESS with long-range order.
However, the efficiency of the dissipative cooling mechanism decreases with the system size N and the correlation length eventually becomes zero in the thermodynamic limit N → ∞.
We provided concrete numerical bounds for the dissipative couplings γ and κ, as well as the system size N in order for the dissipative process to perform robustly.
The following picture appears with regard to experiments: It seems that the finite system size N is crucial for the dissipative cooling protocol to generate macroscopic order in the presence of a nonvanishing coupling to a thermal bath. That is, our results indicate that the driving with jump operators Q xy can be competitive only for small N, where the necessary ratio of dissipative couplings γ/κ to achieve long-range correlations is small. These results are consistent with a recent experimental realization of the dissipative cooling protocol using up to N = 4 particles [27] and it is reasonable to expect that similar proposals for dissipative state generation might face the same limiting constraints with respect to system size.
It is a natural question whether our discussion of closed hierarchies for s = 1/2 quantum spin systems can be generalized to arbitrary spin representations. This would provide a unique means of studying the classical limit (s → ∞) of driven open quantum spin systems and possibly other types of dissipative dynamics with distinct properties of the final state.
Additional information on the asymptotic dynamics for quantum dissipative processes can be obtained by investigating the linear response [60] in the vicinity of the nonequilibrium steady state. This would allow us to inquire about the presence of generalized fluctuationdissipation relations (see, e.g., Refs. [61][62][63]) in the general setting of open Markovian quantum dynamics. We hope to address these questions in future work.