Full decoherence induced by local fields in open spin chains with strong boundary couplings

We investigate an open $XYZ$ spin $1/2$ chain driven out of equilibrium by boundary reservoirs targeting different spin orientations, aligned along the principal axes of anisotropy. We show that by tuning local magnetic fields, applied to spins at sites near the boundaries, one can change any nonequilibrium steady state to a fully uncorrelated Gibbsian state at infinite temperature. This phenomenon occurs for strong boundary coupling and on a critical manifold in the space of the fields amplitudes. The structure of this manifold depends on the anisotropy degree of the model and on the parity of the chain size.


I. INTRODUCTION
Manipulating a quantum system in non-equilibrium conditions appears nowadays one of the most promising perspectives for proceeding our exploration of the intrinsic richness of quantum physics and for obtaining an insight on its potential applications [1][2][3]. In particular, much attention has been devoted to the study of the nonequilibrium steady state (NESS) in quantum spin chains, coupled to an environment, or a measuring apparatus. This is described, under Markovianity assumptions [4][5][6], in the framework of a Lindblad Master equation (LME) for a reduced density matrix, where a unitary evolution, described via Hamiltonian dynamics, is competing with a Lindblad dissipative action. Under these conditions, quantum spin chains subject to a gradient evolve towards a NESS, where spin and energy currents set in. In quasi onedimensional systems, such currents exhibit quite exceptional properties like scalings, ballisticity and integrability [7][8][9][10][11][12][13]. Many of these unexpected features stem from the fact that the NESS, corresponding to a fixed point of the LME dissipative dynamics with a gradient applied at the chain boundaries, are not standard Gibbs-states. Moreover, further peculiar regimes appear when the time lapse between two successive interactions of the quantum chain with the Lindblad reservoir becomes infinitely small, while the interaction amplitude is properly rescaled. In the framework of projective measurements, this kind of experimental protocol corresponds to the so-called Zeno effect, that determines how frequent projective measurements on a quantum system have to be performed in order to freeze it in a given state [14,15].
In this paper we shall rather focus on a Zeno regime for non-projective measurements, that has been found to describe new counterintuitive scenarios for NESS. In particular, in [17] it was shown that in a boundary driven XXZ spin chain, for suitable values of the spin anisotropy the NESS is a pure state. We want to point out the importance of this result in the perspective of engineering dark states, that have the advantage to be more stable against decoherence, than isolated quantum systems and, therefore, better candidates for technological ap-plications [3,16]. Here we investigate how this non projective Zeno regime can be manipulated by the action of a strictly local magnetic field, whose strength is of the order of the exchange interaction energy of the XYZ Heisenberg spin chain model. The main result of our investigations is that, by such a local effect, one can kill any coherence of the NESS and turn it into a mixed state at infinite temperature. More generally, the von Neumann entropy of the NESS can be changed from its minimum value to its maximum one just by tuning the local magnetic field, provided the coupling with the baths is sufficiently strong.
The paper is organized as follows. In Section II we describe the general properties of the non projective Zeno setup and the way the spin XYZ chain is coupled to the Lindblad reservoirs. The effect of complete decoherence induced by the addition of a fine-tuned local magnetic field acting on the spins close to the boundaries are discussed in Section III. A short account of the symmetries characterizing the NESS in the special case of a XXZ spin chain is reported in Section V. In Section VI we investigate the non-commutativity of the different limits to be performed in the model and the presence of corresponding hierarchical singularities. We conclude with a discussion on the perspectives of our investigations (see Section VII).
Appendices A,B, C and D contain some relevant technical aspects.

II. THE MODEL
We study an open chain of N quantum spins, represented by the Hamiltonian operator H , in contact with boundary reservoirs. The time evolution of the reduced density matrix ρ is described by a quantum Master equation in the Lindblad form [4][5][6] (we set = 1) where L L [ρ], L R [ρ] are Lindblad dissipators acting on spins at the left and right boundaries of the chain, respectively. This is an usual setup for studying transport in quantum spatially extended systems, where the explicit choice of L L and L R is suggested by the kind of application one has in mind. In this way, one describes an effective coupling of the chain, or a part of it, with baths or environments. Within the quantum protocol of repeated interactions, Eq.(1) describes an exact time evolution of the extended quantum system, provided the coupling with the Lindblad reservoirs is suitably rescaled [6].
Here we are interested to explore the strong coupling condition, i.e. Γ → ∞, that corresponds to the so-called Zeno regime. In this case one can obtain the stationary solution of Eq.(1) in the form of the perturbative expansion where ρ N ESS (ξ, Γ) is the density matrix of the non equilibrium steady-state and ξ is a symbol epitomizing the model parameters (e.g. bulk anisotropy, exchange energy, magnetic field, etc.). Suppose that the stationary solution ρ N ESS (ξ, Γ) is unique. This fact will be validated further for all our examples. Moreover, the first term of expansion (2), i.e. ρ 0 = lim Γ→∞ lim t→∞ ρ(Γ, ξ, t), satisfies the stationarity condition L LR [ρ 0 ] = 0, where L LR = L L + L R is the sum of the Lindblad actions in (1) . This suggests that ρ 0 can be represented in a factorized form where ρ L and ρ R are the one-site density matrices at the chain boundaries, satisfying L L [ρ L ] = 0 and L R [ρ R ] = 0, and M 0 is a matrix to be determined self-consistently. It is convenient to separate explicitly the identity matrix I 2 ⊗N−2 from M 0 , in such a way that M 0 is a traceless operator, due to the condition T r(ρ 0 ) = 1 . By substituting the perturbative expansion (2) into Eq.(1) and by equating terms of the order Γ −k , one can easily obtain the recurrence relation whose general solution has the form 5) provided the following secular conditions (for more details see [25]) are satisfied where ker(L LR ) denotes the nucleus of the operator L LR . Notice that, in order to obtain an explicit solution, one has to compute the inverse operator L −1 LR , that appears in Eq.(5). In summary, Eqs (3), (5) and (6) define a general perturbative approach, that applies in the Zeno (i.e., strong coupling) regime.

We consider the Hamiltonian
is the Hamiltonian of an open XY Z Heisenberg spin chain and V l is a local inhomogeneity field acting on spin l to be specified later on (see . Moreover, we consider Lindblad dissipators, L L and L R , favouring a relaxation of boundary spins at k = 1 and k = N towards states described by one-site density matrices ρ L and ρ R , i.e. L L [ρ L ] = 0 and L R [ρ R ] = 0. In particular, we choose boundary reservoirs that tend to align the spins at the left and right edges along the directions l L and l R , respectively. These directions are identified by the longitudinal and azimuthal coordinates as follows: Such a setting is achieved by choosing the Lindblad action in the form L and In the absence of the unitary term in (1), the boundary spins relax with a characteristic time Γ −1 to specific states described via the one-site density matrices The condition L A [ρ A ] = 0 follows from definition (8), while the relations S A S † A = ρ R and (S A ) 2 = (S † A ) 2 = 0 can be easily checked.
In analogy with [18], it can be easily shown that, for the chosen boundary dissipation setup described by Eqs (8), (9) and (10), the NESS is unique. By applying the perturbative approach in the Zeno regime, one finds that the unknown matrices M k (∆) are fully determined by secular conditions (6). As shown in Appendix A, for the specific choice (8) of the Lindblad operators, they are equivalent to the requirement of a null partial trace We want to point out that the computation of the full set of matrices {M k (∆)} for any ∆ = 0 is quite a nontrivial task. However, in the Zeno limit, Γ → ∞, we are just interested in computing the zero-th and the first order contributions M 0 , M 1 , which can be completely determined by solving the set of secular equations (13) for k = 0, 1, 2.

III. MANIPULATIONS OF NESS BY NON UNIFORM EXTERNAL FIELDS
The properties of the model introduced in the previous section have been widely investigated for V l = 0 and ϕ = π/2 in [25]. Here we are interested in studying how the properties of the NESS can be modified when V l is an additional local field, that corrupts the homogeneity of the XYZ spin chain.
Notice first that a local field applied to the boundary spins at positions k = 1 and k = N does not affect the strong coupling limit ρ 0 = lim Γ→∞ ρ N ESS (Γ). On the other hand, applying a local field to the spins at positions k = 2 and k = N − 1 can modify ρ 0 in a nontrivial way. The Hamiltonian reads where Carrying out the procedure outlined in the previous section, we can find the form of the density matrix of the NESS in the Zeno regime, ρ 0 . This is a function of the angles θ L , ϕ L , θ R , ϕ R , of the anisotropy parameter ∆ and of the local fields h, g. One can argue that, in general, the NESS should be an entangled state, depending in a nontrivial way on the local fields. Due to the boundary drive, the NESS typically exhibits nonzero currents (magnetization current, heat current, etc.), irrespectively of the presence of the local fields. However, in the Zeno limit, there are critical values of the local fields for which a complete decoherence of the NESS occurs.
More precisely, we formulate our results under the following boundary condition assumptions: • targeted boundary polarizations are neither collinear nor anti-collinear ( l L = ± l R ); • at least one of the polarizations (e.g. the left targeted polarization) is directed along one of the anisotropy axis X, Y, or Z; • the corresponding local fields ( h at site 2 and g at site N − 1) are collinear to the respective targeted boundary Then, there exists a zero-dimensional or a one dimensional critical manifold in the h, g-plane (h cr , g cr ), such that, in the Zeno limit , the NESS on this manifold becomes Notice that this a peculiar state: apart the frozen boundary spins, all the internal spins are at infinite temperature. Indeed, tracing out the boundary spins, one obtains the Gibbs state at infinite temperature Also notice that on the critical manifold the Von-Neumann entropy of the NESS, S V N E = −T r(ρ N ESS log 2 ρ N ESS ), in the Zeno limit attains its maximum value given by since ρ L , ρ R are pure states. In the following, we also refer to state (17) as the state of maximal decoherence. We have performed explicit calculations (see below) that confirm the above statement for different spin chains up to N = 8. The particular form of the NESS assumed in these cases, however, strongly suggests that the above results maybe of general validity and the critical manifold (h cr , g cr ) independent on N .
The critical manifold has been fully identified for the following cases.
-XYZ chain: J x = J y = ∆. If the left, l L , and the right, l R , polarizations point in directions of different principal axes where e X = (1, 0, 0), e Y = (0, 1, 0), e Z = (0, 0, 1), then for chains with an even number, N , of spins, the manifold (h cr , g cr ) consists of three critical points: For odd N , the critical point P α,β is missing and the critical manifold reduces only to the points P α , P β , above. If only one of the two boundary driving points in the direction of a principal axis, the critical manifold reduces to a single point, either P α or P β , for both even and odd N .
-XXZ chain: J x = J y = J = ∆. If both l L and l R lay onto the XY -plane, we can parametrize the targeted boundary polarizations via a twisting angle in the XY -plane ϕ as θ 1 = θ 2 = π 2 , ϕ 1 = ϕ, ϕ 2 = 0, corresponding to l L = (cos φ, sin φ, 0) and l R = (1, 0, 0). The critical fields are aligned parallel to the targeted boundary magnetization, i.e. h cr = (h cr cos φ, h cr sin φ, 0) , g cr = (g cr , 0, 0), and we find the one-dimensional critical manifold Notice that this expression is independent of system size N , of the anisotropy ∆ and of the twisting angle ϕ. If one of the two targeted polarizations points out of the XY -plane, the critical manifold becomes zero-dimensional and consists of one, two or three critical points (depending on the polarization direction and on N being even or odd) as discussed for the full anisotropic case.
-XXX chain: J x = J y = ∆ ≡ J. The critical manifold for arbitrary non-collinear boundary drivings is one-dimensional and it is given by Eq. (20).
The above statements are illustrated in Figs. 1 and 2 for the case of a chain of N=4 spins. In particular, in Fig. 1 we show a contour plot of the VNE surface as a function of the applied fields for the XYZ case with left and right boundary polarizations fixed along the X and Z directions, respectively. where the VNE reaches its maximum value SV NE = 2 and the NESS becomes a completely mixed state, respectively. Other parameters are fixed as lL = eX, lR = eZ. Green, yellow, pink, orange, brown, red and blue contour lines refer to SV NE values: 1, 1.2, 1.4, 1.6, 1.8, 1.9, 1.95, respectively. Notice the presence of the narrow corridors around PX and PZ in which the deviation, SV NE −2, of the VNE from its maximum value becomes very small.
The three critical points P X , P Z , P X,Z mentioned above correspond to the green, red, and white dots shown in the top panel of the Figure. Notice the presence of narrow corridors (blue shaded) around the P X and P Z critical points, inside which the VNE keeps very close to the maximal value S V N E = 2 but never reach it, except at the critical points. This is quite different from the partially anisotropic XXZ case shown in Fig. 2, where the existence of the critical line (blue line) is quite evident.
Similar results are found also for longer chains. In particular, in Fig.3 we show a cut of the VNE surface for a partially anisotropic XXZ chain of N = 5 spins. For the sake of simplicity we have set J x = J y = 1 and considered the cut at h = 0 so that the VNE of the NESS, in the Zeno limit, becomes a function of g only. We see that for g = g cr = −2, the VNE reaches the maximum value N − 2 indicating that the corresponding NESS has the form (17).
As to the dependence of the critical manifold on parity of N , we find that while for odd sizes N = 3, 5, 7 and XY Z Hamiltonian ( see Fig. 4, top panel for an illustration) there are only two critical points (the critical point P α,β is missing), for even sizes N = 4, 6, 8 cases there are three critical points. These observations strongly suggest a qualitative difference between even and odd N in the model, which is manifested in other NESS properties as well, see e.g. (22), (23).
It is worth to note here that for h = g = −J, i.e. the case excluded in (20), the NESS behaves non-analytically in the Zeno limit Γ → ∞. As we are going to discuss in Sec.VI, this non-analyticity is a consequence of the non-commutativity of the limits Γ → ∞ and h = g → −J.
Conversely, for any finite boundary coupling Γ, i.e. far from the Zeno limit, the NESS is analytic for arbitrary amplitudes of the local fields (the first order correction to the NESS for large Γ is proportional to Γ −1 as shown in the App.C). This is also seen from Fig.5 where the VNE of the NESS is reported as a function of the local field g for different values of the boundary coupling Γ and same parameters as in Fig.3 (see curve ∆ = 0.6). Notice that the thin black line obtained for Γ = 10 3 , is already in full overlap with the Zeno limiting curve depicted in Fig.3 for ∆ = 0.6. Also note the persistence of the peak at g = −2 even for relatively small values of Γ away from the Zeno limit. Similar behaviors are observed for different choices of boundary polarizations and of local fields (not shown for brevity), thus opening the possibility to detect the signature of the above phenomena in real experiments. In this respect we remark that the near-boundary magnetic field h and the anisotropy ∆ as suitable parameters for controlling the dissipative state of the system in a NESS. Thus, if g = 0, the NESS can be made a pure state by tuning the anisotropy ∆ to a specific value ∆ * (ϕ, N ). For instance, we find that for g = 0 and ∆ * ± (π/2, 5) = 1 2 ± 1 2 √ 2 the NESS is a pure state [26], while for g cr = −2, the NESS in the bulk becomes an infinite temperature state (17), i.e. a maximally mixed state. Thus, by suitably tuning the anisotropy and the local magnetic field one can pass from minimally mixed (pure) NESS state to a maximally mixed one.
It should be emphasized at this point that general thermodynamic equilibrium quantities, e.g. the temperature, are not well-defined for a generic NESS. In fact, pure states allowed by Liouvillian dynamics are not ground states of the Hamiltonian, but are characterized by a property of being common eigenvectors of a modified Hamiltonian and of all Lindblad  operators [16,21]. Likewise, an absence of currents in the NESS does not necessarily imply a thermalization of the system: in fact also for fully matching boundary conditions the NESS is not a Gibbs state at some temperature, so that correlation functions remain far from those of an equilibrium system. From this point of view, the decoherence effect described in present paper can be viewed as a reaction of a nonequilibrium system on a local perturbation (the local magnetic field): as is well-known, a local perturbation in nonequilibrium can lead to global changes of a steady state.
On the other hand, a fully mixed state as such has appeared already in the context of driven spin chains: if both Lindblad boundary reservoirs target trivial states with zero polarization (ρ R = ρ L = I/2), the NESS is maximally mixed ρ N ESS = I 2 ⊗N , which is a trivial solution of the steady Lindblad equation for any value of boundary coupling. The respective NESS is often being referred to as a state with infinite temperature [22]. Note that our case is drastically different from the latter: the maximally mixed state (17) appears only in the bulk, after tracing the boundary spins, in a system with generically strong boundary gradients, and under strong boundary coupling.
A few more remarks are in order: (i) the amplitudes of the critical local fields scale with the amplitude of the Hamiltonian exchange interaction, i.e. h cr → γh cr if H XXZ → γH XXZ (this is a consequence of the linearity of the recurrence relations (5) and (13) ); (ii) the NESS may take the form (17) only in the Zeno limit Γ → ∞; in fact, the first order correction to the NESS is proportional to Γ −1 and does not vanish (see Appendix C); the fully decoherent state (17) is intrinsic to nonequilibrium conditions and, strikingly enough, it persists even for nearly matching or fully matching boundary driving, as we are going to discuss in Sec.IV.
We want to conclude this Section by pointing out that a fully analytic treatment of the problem for arbitrary large values of N should encounter serious technical difficulties. The main one concerns the solution of the consistency relations determined by the secular conditions (6) for the perturbative expansion (2), with the zero-order term given by (17). Moreover, finding the first order correction to NESS, proportional to Γ −1 , amounts to solve a system of equations, whose number grows exponentially with N . With Matematica we were able to solve that system of equations analytically for N ≤ 5 and numerically up to N ≤ 7.

IV. MATCHING AND QUASI-MATCHING BOUNDARY DRIVINGS
In the previous Section we have discussed the case where the complete alignment of the boundary Lindblad baths was excluded. In this Section we want to analyze the specific case where they are aligned (or quasi-aligned) in the same direction on the XY -plane.
A complete alignment, i.e. l L = l R , corresponds to a perfect matching between the left and right boundary Lindblad baths, that yields a total absence of boundary gradients, so that any current of the NESS vanishes. Also in this case the Gibbs state at infinite temperature can be achieved by suitably tuning the values of the near-boundary fields, but for even-sized chains, only.
Let us first illustrate this finding for the XYZ case. With no loss of generality, we can set l L = l R = e Z = (0, 0, 1). The behavior of the driven chain with local fields depends drastically on whether the size of the chain N is an even or an odd number: in the former case we find the critical one dimensional manifold, defined by h cr + g cr = −2∆, h cr = −∆ ; (21) in the latter case N = 3, 5, .., we do not find any critical point. This result has been found explicitly for 3 ≤ N ≤ 6, but, since it depends on the effect of local perturbations, it seems reasonable to conjecture that it should hold for larger finite values of N . This result holds as long as the Heisenberg exchange interaction in the plane perpendicular to the targeted direction (the XY -plane in this example) is anisotropic, i.e. J x = J y . Conversely, for J x = J y , the infinite temperature state (17) cannot be reached for any value of the local fields h and g. There is a delicate point to be taken into account when we fix h = h cr and we perform the limit J y → J x , i.e. we reestablish the model isotropy: for complete alignment, l L = l R = e Z , the NESS is singular. The way this singularity sets in is shown in Fig. 6. In the limit when the anisotropy in the direction transversal to the targeted direction becomes infinitesimally small |J y − J x | → 0 the NESS is a pure state with minimal possible S V N E → 0 for any amplitude of the local field values, except at a critical point where S V N E is maximal. The noncommutativity of similar limits and the dependence of the NESS properties on the parity of system size N in Lindblad-driven Heisenberg chains, have been observed in previous studies [23], [24]. Also in these cases, the origin of noncommutativity is a consequence of global symmetries of the NESS, that, for our model, are discussed in Section V. In the isotropic case, as long as the local fields are parallel to the targeted spin polarization, the NESS does not depend on them: it is a trivial factorized homogeneous state with a maximal polarization matching the boundaries, i.e. ρ N ESS = (ρ L ) ⊗N . This can be easily verified by a straightforward calculation.
Another kind of NESS singularity can be found in the partially anisotropic case, with quasi-matching boundary driving in the XY isotropy plane. As a mismatch parameter we inroduce the angular difference between the targeted polarizations at the left and the right boundaries, ϕ = ϕ L − ϕ R . For ϕ = 0 and in the absence of local fields, we have found that the spin polarization at each site of the chain is parallel to the targeted polarization; on the other hand, even in the Zeno limit, it does not saturate to the value imposed at the boundaries j = 1 , N . In general, this is not an equilibrium Gibbs state, even in the Zeno limit and for any finite boundary coupling Γ. However, if the near-boundary fields are switched on and tuned to their critical values, the coherence of this state is destroyed and the NESS becomes an infinite temperature Gibbs state. On the other hand, we have found that there is a relevant difference between quasi-matching and mismatching conditions for even and odd values of N (notice that the isotropic and the free fermion cases, ∆ = 1, ∆ = 0, are special and should be considered separately). Our results can be summarized as follows: -N odd. We can fix the boundary mismatch by choosing ϕ L = ϕ, ϕ R = 0, the left local field h = 0, and study the NESS as a function of the right local field g. At g = g cr = −2, the NESS becomes trivial (maximally mixed); however, as shown in panel (a) of Fig.7, for small mismatch we find a singular behavior of the NESS close to g = g cr . Analytic calculations (not reported here) show that for ϕ = 0 there is a singularity at g = g cr , as a result of the non-commutativity of the limits ϕ → 0 and g → g cr .
-N even. Unlike the previous case, the NESS is analytic for  Fig.7). For g = g cr the NESS becomes trivial (maximally mixed), also for ϕ = 0. Finally, let us comment about two special cases, for "equilibrium" boundary driving conditions, i.e. ϕ L = ϕ R . For ∆ = 0 (free fermion case), the NESS is a fully mixed state (apart from the boundaries) for all values of g. For ∆ = 1 (isotropic Heisenberg Hamiltonian), the NESS is a trivial factorized state, fully polarized along the axis of the boundary driving, for any value of g. Both statements can be straightforwardly verified.
NESS singularities, onset of which can be recognized in Fig. 6 and Fig. 7a, appear because of non-commutativity of limits. Noncommutativity of various limits, implying singularities and nonergodicity, which are due to global symmetries is a well-established phenomenon and occurs already in Kubo linear response theory describing fluctuations of a thermalized background. In nonequilibrium open quantum systems, however, the presence of NESS symmetries at special value of parameters is manifested much strongly, due to richer phase space which includes both bulk parameters (such as anisotropy and external field amplitudes) and boundary pa-rameters (such as coupling strength). As a result, noncommutativity of the limits and consequent NESS singularities seems to be a rather common NESS feature. In the next two sections we reveal some of NESS symmetries and show that the respective singularities, connected with them, can be observed already in a finite system consisting of a few qubits.

V. SYMMETRIES OF NESS
Symmetries of the LME are powerful tools that reveal general, system size-independent properties of the Liouvillean dynamics (1). In the case of multiple steady states, symmetry based analysis allows one to predict different basins of attraction of the density matrix for different initial conditions [19]. For a unique steady state, symmetry analysis provides a qualitative description of the Liouvillean spectrum [20] or the formulation of selection rules for steady state spin and heat currents [23]. It is instructive to list several general NESS symmetries valid for our setup. We restrict to XXZ Hamiltonian with J x = J y = 1, and perpendicular targeted polarizations in the XY -plane, i.e. l L = (0, −1, 0), l R = (1, 0, 0). The LME has a symmetry, depending on parity of N , which connects the NESS for positive and negative ∆. Let us denote by ρ N ESS (N, ∆, h, g, Γ) the nonequilibrium steady state solution of the Lindblad master equation (see (1) and (8) ) for the Hamiltonian (B1) reported in Appendix B. It is known that this NESS is unique [18] for any set of its parameters; moreover, one can easily check that ⊗ σ z n and the asterisk on the r.h.s. of both equations denotes complex conjugation in the basis where σ z is diagonal. Eqs (22) and (23) hold for any value of the local fields h, g and for any coupling Γ, including the Zeno limit Γ → ∞. Due to properties (22) and (23), we can restrict to the case ∆ ≥ 0 further on. For g = −h, ρ N ESS (N, ∆, h, g, Γ) has the automorphic symmetry, where R(A⊗B ⊗...⊗C) = (C ⊗....⊗B ⊗A)R is a left-right reflection, U rot = diag(1, i) ⊗N is a rotation in XY plane, U rot σ x n U + rot = σ y n , U rot σ y n U + rot = −σ x n , and Σ x = (σ x ) ⊗N .
Here we consider the XXZ Hamiltonian and a perpendicular targeted polarizations in the XY -plane l L = (0, −1, 0), l R = (1, 0, 0); the near-boundary fields are taken on the critical manifold, i.e. h + g = −2. For N = 3, 5 and ∆ > 0 we have found the noncommutativity conditon Making use of (17), the r.h.s. of (25) can be rewritten (26) For the simplest nontrivial case N = 3, the validity of these noncommutativity relations is verified by the calculations reported in Appendix B (see (B4)). On top of (25), we find an additional singularity at the isotropic point ∆ = 1 for N > 3 Due to the symmetry conditions (22) and (23), the singularity is present also for ∆ = −1. Eqs (25) and (27) , engendering a zero-dimensional submanifold of the critical manifold. Thus, a hierarchy of singularities is formed. It is quite remarkable that such hierarchical singularities emerge without performing the thermodynamic limit N → ∞. In fact, as shown in Appendix D, they can be explicitly detected already for N = 4. For N = 5 we have found other singular manifolds, parametrized by h, g, and ∆. For the sake of space, details of this case will be reported in a future publication. The appearance of the singularity at h → −1, g → −1 is a consequence of the additional symmetry (24) at this point. By direct inspection of the analytic formulae obtained for N = 3, 4, 5, we can guess the form of the limit state lim Γ→∞ lim h→1 lim ∆→1 as a fully factorized one, namely Conversely, for generic ∆ and odd N ≥ 5, the limit state lim Γ→∞ lim h→−1 ρ N ESS (N, h, ∆, Γ) does not take a factorized form. Notice also that from making use of Eqs (22) and (23), we readily obtain also the NESS limit state for ∆ → −1:

VII. CONCLUSION.
In this paper we extensively analyzed the properties of the NESS of open Heisenberg spin chains, subject to the action of LME at their boundaries and of perturbing magnetic fields at the near-boundary sites. The setup we deal with operates in the Zeno regime, i.e. in the strong coupling limit, Γ → +∞ (see Eq. (1) ). Most of our analytic and numeric calculations have been performed for relatively small values of the chain size N . On the other hand, as a consequence of the local nature of the reservoirs and of the perturbing magnetic fields, we conjecture that many of these results could be extended to large finite values of N : the delicate question of how they might be modified in the thermodynamic limit is still open. At the present level of standard computational power, the strategy of performing large scale calculations to get any inference on such a limit is impractical, because the number of equations to be solved grows exponentially with N .
Despite all of these limitations, the main outcome of our study is quite unexpected: by tuning the near-boundary magnetic fields we can manipulate the NESS, making it pass from a dark pure state (for a suitable choice of the value of the anisotropy parameter ∆), to a fully uncorrelated mixed state at infinite temperature.
We have also discussed how this general scenario emerges in the anisotropic, partially anisotropic and isotropic cases. The influence of different alignment conditions imposed by the Lindblad reservoirs has been extensively explored, together with the symmetries of the NESS and their importance for engendering hierarchical singularities due to the noncommutativity of different limits, performed on the model parameters.
A physically relevant point in our discussion concerns the possibility of performing such a manipulation of the NESS also for large but finite values of Γ: numerical investigations confirm this expectation, thus opening interesting perspectives of experimental investigations.
Noting that Q has the property 1 2 L LR Q = −Q, we obtain from (4) and (5) the first order correction to ρ 0 Let us assume that M 1 = 0. Then, in the second order of perturbation theory, we have After some calculations we obtain where the unwanted secular terms are written out explicitly, and T r 1,N R = 0. The unwanted terms proportional to ∆ do not depend on h, g. For any ∆ = 0 the secular conditions T r 1,N [H, ρ 1 ] = 0 cannot be satisfied. This contradiction shows that M 1 = 0 for any ∆ = 0.