Nonequilibrium Quantum Phase Transitions in the XY model: comparison of unitary time evolution and reduced density matrix approaches

We study nonequilibrium quantum phase transitions in XY spin 1/2 chain using the $C^*$ algebra. We show that the well-known quantum phase transition at magnetic field $h = 1$ persists also in the nonequilibrium setting as long as one of the reservoirs is set to absolute zero temperature. In addition, we find nonequilibrium phase transitions associated to imaginary part of the correlation matrix for any two different temperatures of the reservoirs at $h = 1$ and $h = h_{\rm c} \equiv|1-\gamma^2|$, where $\gamma$ is the anisotropy and $h$ the magnetic field strength. In particular, two nonequilibrium quantum phase transitions coexist at $h=1$. In addition we also study the quantum mutual information in all regimes and find a logarithmic correction of the area law in the nonequilibrium steady state independent of the system parameters. We use these nonequilibrium phase transitions to test the utility of two models of reduced density operator, namely Lindblad mesoreservoir and modified Redfield equation. We show that the nonequilibrium quantum phase transition at $h = 1$ related to the divergence of magnetic susceptibility is recovered in the mesoreservoir approach, whereas it is not recovered using the Redfield master equation formalism. However none of the reduced density operator approaches could recover all the transitions observed by the $C^*$ algebra. We also study thermalization properties of the mesoreservoir approach.


Introduction
Equilibrium phase transitions are determined as non-analyticities of free energy and can appear only in the thermodynamic limit [1]. At finite temperatures, phase transitions are driven by thermal noise, and states of systems are minima of the free energy; in other words, systems tend to show large entropy for higher temperatures. On the other hand, quantum phase transitions are driven by quantum fluctuations and appear at absolute zero temperature [2]. While the entropy contribution to the free energy in equilibrium systems prohibits the equilibrium phase transition in one-dimensional systems with short-range interactions [3,4], a nonequilibrium phase transition in one dimension is possible. Nonequilibrium phase transitions are usually considered as qualitative changes of the steady state. In classical cases, nonequilibrium phase transitions are well studied and appear, for example, in driven diffusive models [5]. The Yang-Lee description of equilibrium phase transition in terms of the zeros of the partition function [1] has been applied to a classical nonequilibrium setting [6].
On the other hand, nonequilibrium quantum phase transitions (NQPT) are much less known. In that respect, some interesting numerical results [7,8] show an NQPT in onedimensional boundary-driven nonequilibrium spin systems. Moreover, approaches to deal with the nonequilibrium quantum systems usually involve different approximations, and their validity near an NQPT has not been carefully discussed. Thus, it is important to analyze the differences among several approaches to nonequilibrium quantum systems in a simple pedagogic model undergoing an NQPT.
Nonequilibrium steady states (NESS) of quantum systems are generally studied using two approaches.
• One approach is to decompose the total, composite system into a finite system and reservoir parts and derive, by tracing out the latter, an effective master equation for the density operator of the former; i.e., we obtain an effective master equation for the reduced density operator of the finite system. It is possible to derive an exact master equation for the reduced density operator of the system, however, this is usually as difficult to solve as the original problem. Therefore, various approximation schemes have been developed to suitably describe the different regimes. In the simplest case, we obtain a completely positive Markovian master equation, namely the Lindblad master equation [9,10], and the NESS is obtained as a projection of the initial density operator onto the null-space of the generator of the dynamics, the Liouvillian  . Markovian master equations have mostly been used in quantum optics [11], quantum information, and quantum computation as a simple model of noisy channels. However, recently they have also been applied in the context of condensed matter to study the high-temperature transport properties of simple one-dimensional systems [12][13][14][15] and to discuss special nonequilibrium states of matter [16]. In this paper, we shall employ two models, which have been developed to study transport of open quantum systems, namely the modified Redfield master equation approach [8] and the mesoreservoir approach [17]. • The other is the * C algebra approach, which was first introduced with the purpose of rigorously formulating equilibrium statistical mechanics [18]. It has been applied to infinitely extended systems [19]. Starting from Ruelleʼs work [20,21] on scatteringtheoretical characterizations of NESS, the * C algebra method has been extensively developed (see, for instance, references of [22]) and also applied in the context of phase transitions [23,24]. Contrary to the reduced density operator method, the * C algebra deals with infinite systems which evolve unitarily. It was shown that if a system is connected to two separated subsystems initially in thermal states, a finite part of the total system approaches a unique NESS under some mathematical conditions [20]. In that case, the constructed NESS does not depend on the initial decomposition, i.e., the size of the central (finite) system, or the initial density operator of the central system. Although the * C algebra provides mathematically rigorous results, the applications are quite limited. Accordingly, it is beneficial to use the * C algebra approach to test other approaches for their validity.
In this paper we use the * C algebra approach, the Lindblad formalism with mesoreservoirs, and the modified Redfield master equation to study the XY spin 1/2 chain-a paradigmatic model exhibiting quantum phase transitions (QPTs) [25][26][27]. The Hamiltonian of the transverse-field XY spin 1/2 chain is where γ denotes the anisotropy, h denotes the magnetic field, and σ m x,y,z are the Pauli spin operators at the mth site. To begin with, let us summarize the phase transitions of the XY spin 1/2 chain discussed previously. At h = 1, there is a second order QPT characterized by the order parameter σ σ 〈 〉 α α l m (α = x for γ > 0 and α = y for γ < 0), which separates the ferromagnetic ordered phase ( < h 1) and the paramagnetic disordered phase ( > h 1). The equilibrium average of σ 〈 〉 α m for α = x, y is always zero, and it cannot be used as an order parameter. Magnetization in the z direction is always finite, but susceptibility is divergent at absolute zero temperature and h = 1. In addition, at γ = 0 and | | < h 1 there is a quantum phase transition between the ordered phase in the x direction (γ > 0) and the y direction (γ < 0). In case of the open XY spin 1/2 chain coupled to local Lindblad reservoirs [7] and Redfield reservoirs [8], indications of NQPT were reported. Numerically, it was observed that γ = − h 1 c 2 is a critical magnetic field which separates phases showing different scaling of the quantum mutual information (QMI), spectral gap, and far-from-diagonal spin-spin correlations σ σ . This transition was recently described by an information-geometric approach using a fidelity distance measure [28].
We shall complement the already-obtained phase diagram of the XY chain by using the * C algebra method. We focus on the NQPT and reveal new quantum phase transitions, which do not exist in equilibrium, and show the existence of an NQPT at h = 1 for arbitrary temperatures of reservoirs, which is associated with a discontinuity of the third derivative of the off-diagonal elements of the correlation matrix. We also demonstrate a discontinuity of the first derivative of the correlation matrix elements at = h h c . In addition, we show that QPT at h = 1 persists in the nonequilibrium case if the temperature of at least one reservoir is set to be zero. A phase diagram showing equilibrium and nonequilibrium quantum phase transitions of the XY spin 1/2 chain is presented in figure 1.
We shall compare NQPTs obtained using the * C algebra method with those obtained by simulations of the reduced density operator approaches, namely the mesoreservoir method [17] and the modified Redfield master equation [8]. We also discuss the thermalization properties of the mesoreservoir approach.
The rest of the paper is structured as follows. In section 2 we apply the * C algebra approach and analytically show the four different nonequilibrium quantum phase transitions of the XY model. In section 3 we discuss the master equation approaches. In subsection 3.1, we first present the results of the Lindblad mesoreservoir approach and discuss its equilibrium and nonequilibrium properties. Next, in subsection 3.2, we study the modified Redfield model. We discuss the results and conclude in section 4.
2. NQPT in the XY spin 1/2 model: the C Ã algebra approach In this section we study the NQPTs of the XY spin 1/2 chain using the results of the * C algebra approach. We shall show the coexistence of two NQPTs at h = 1 and γ ≠ 0 (the critical XY line) and an NQPT at the critical magnetic field γ = ≡ − h h 1 Figure 1. Phase diagram of the equilibrium and nonequilibrium XY spin 1/2 model. Quantum phase transitions: at γ = | | < h 0, 1 (critical XX) there is a transition from the ordered phase in the x direction for γ > 0 and the ordered phase in the y direction for γ < 0. At h = 1 (critical XY) the system has divergent magnetic susceptibility in the z direction. Nonequilibrium quantum phase transitions: at the nonequilibrium critical line there is a jump in ∂ 〈 〉 The XY spin 1/2 model has been extensively studied with the aid of the * C algebra; however, the nonequilibrium quantum phase transitions have not been discussed so far. We briefly explain the terminology of the * C algebra method in appendix B. In [29,30] the NESS of the model was rigorously constructed using the scattering theory proposed by Ruelle [20]. Long-range correlation [29] and non-negativity of entropy production [21] have been discussed for this NESS. Those works are highly mathematically oriented. In contrast, we provide a simpler and more direct calculation, and focus on the nonequilibrium phase transitions. We start with an infinite chain, which is separated into the left semi-infinite part (L), the central finite part (S), and the right semiinfinite part (R). The infinite chain is initially in the product state ρ where ρ S is an arbitrary state and ρ L,R is the density operator of the canonical ensembles with temperature T L,R (ρ ∼ − e H T / L,R L,R ). One can prove that a unique NESS of the chain exists and that the diagonal modes with positive (negative) velocities are distributed according to the canonical ensemble of the left (right) reservoir, and satisfy Wickʼs theorem, which is explicitly expressed as equations (5) and (6).
Let us first diagonalize the XY spin 1/2 model. With the aid of the Araki-Jordan-Wigner transformation (see appendix B), the Hamiltonian is mapped to a chain with fermions The Hamiltonian (3) can easily be diagonalized  The NESS of this model uniquely exists and is fully characterized by [29,30] and Wickʼs theorem [20], where 〈 · 〉 represents a NESS average, ϵ = + To be more precise, modes with positive and negative velocity follow different KMS conditions. (See equation (27) of [31] for the XY chain. See also equation (54) in [30] for the XX chain, where a more explicit expression in terms of the Fermi distribution is provided in equation (65).) Let us briefly discuss the intuitive idea behind equations (5) and (6), which is implemented rigorously with the * C algebra approach. Note that in the diagonal form (4) the Hamiltonian can be interpreted as a sum of Hamiltonians . Consider a free system such as one described by h 1 . At t = 0 it is split into three parts: a left semi-infinite chain, a central finite chain, and a right semi-infinite chain. The left (right) semi-infinite chain is in a thermal state with temperature T L (T R ). At = + t 0 the three pieces are connected and the particles (or quantum waves described by the diagonal modes) that were confined to the left can now propagate to the right without any scattering; the same holds true for the particles at the right. In that way, every right-going mode in the system ( > v 0 k ) comes from the left and is populated accordingly with ϵ f ( ) k L , and every left-going mode in the system ( < v 0 k ) comes from the right and is populated by . This is the content of equation (5). Equation (6)   Let us first discuss the magnetization that indicates a second-order quantum phase transition (for instance, see the proof of theorem 2.3. in [29]) is the susceptibility in NESS. We shall now show that this transition also persists in the nonequilibrium setting where one of the reservoirs is initiated at finite nonzero temperature.
From equation (8), one can see that the magnetization σ 〈 〉 m z does not depend on the spatial variable m, and both the magnetization and susceptibility are simply the average of those in equilibrium. Due to , the susceptibility diverges if and only if at least one of the reservoirs has absolute zero temperature. The divergence with respect to temperature ν T and the difference of the magnetic field from one | − | h 1 is logarithmic. Therefore, at h = 1 we have a quantum phase transition and an NQPT associated with the divergence of the susceptibility.
The mechanism of the NQPT associated with the logarithmic divergence of the magnetic susceptibility is the same as in the equilibrium case, since the NESS average of the magnetization is a sum of terms coming from the left and right reservoirs in equilibrium. Since all real parts of the correlation matrix have this property, a genuine NQPT should be discussed through 〈 〉 † f f Im l m , which vanishes for the equilibrium state (the kernel is proportional to the difference in Fermi distributions calculated with the initial temperatures of the left and right semi-infinite parts). The first derivative of The presence of the step function θ v ( ) k in (9) produces a discontinuity at = h h c . Physically it can be explained by a change in the direction (sign of velocity) of the diagonal modes c k as a function of the magnetic field h. In the equilibrium case, diagonal modes c k follow the same Fermi distribution independent of the direction (velocity) of the modes. On the other hand, in NESS, diagonal modes c k with positive and negative velocities follow different Fermi distributions. Therefore, changing the direction of the modes alters their nature only in nonequilibrium and this is the origin of the NQPT at = h h c , which is consequently a genuine nonequilibrium phenomenon.
Interestingly, there is also a discontinuity in the third derivative at h = 1. The third From equation (10), one can evaluate the jump of the third derivative at h = 1: Since this transition is not related to the change in the direction of the modes, its origin is different from the transition discussed for = h h c . The behavior of derivatives (9) and (10) is depicted in figure 2. It is interesting that the non-continuities of the derivatives are also present in the case of non-zero temperatures of the initial states of reservoirs ≠ > T T 0 L R , while the transition of a real part such as magnetization is possible if at least one of the reservoirs is set to absolute zero temperature. Thus, the behavior of the imaginary part is quite different from the real part, which is also the reason we believe that 〈 〉 † f f Im l m determines genuine nonequilibrium properties. We recall that imaginary parts are related to currents of conserved quantities, e.g. the energy current and the spin current (the latter is well-defined only in the isotropic case). At γ = 0 and = = h h 1 c we find a square root divergence of magnetic susceptibility with the temperature T L,R and the difference of In this case γ = = h ( 0, 1), the finite temperature nonequilibrium transitions associated with discontinuities in the derivatives of the imaginary parts of the correlation functions disappear. An equilibrium and nonequilibrium phase diagram of the XY spin 1/2 model is depicted in figure 1.
Finally, we discuss scaling of the quantum mutual information (QMI) in the NESS of the * C algebra. The scaling of QMI has recently been discussed in the context of area laws. In [32] it was shown that Gaussian thermal states obey the area law for QMI. On the other hand, the steady states may violate the area law despite being Gaussian. In fact, QMI was used to characterize the NQPT at = h h c in the XY model with Redfield reservoirs [8] (see section 3.2). In [8] it was shown that below the critical field ( = h h c ) QMI obeys the area law, whereas above the critical field QMI scales linearly with the system size. We calculate the scaling of the QMI in the NESS by numerically diagonalizing the two-point Majorana correlation matrix of the * C algebra. The Von Neumann entropy of the block of size n can be expressed in terms of the eigenvalues λ k of the × n n 2 2 two-point Majorana correlation matrix (defined in (14)) as [33] ∑ The quantum mutual information of the block of size n 2 is then given by = − I n S n S n (2 ) 2 ( ) (2 ). In figure 3 we show that the QMI scales logarithmically with the block size n. Interestingly, this scaling does not depend on the system parameters and therefore does not capture the nonequilibrium phase transitions. A related analytical calculation of the area law violation in the * C algebra NESS was recently reported in the free fermion case [34], which is equivalent to the isotropic limit of the XY spin 1/2 chain.

NQPT in the XY spin 1/2 model: master equation approach
In this section we compare the obtained, exact * C algebra results with two models of open quantum systems, namely the Lindblad mesoreservoir and the modified Redfield master equation. In both cases, the dynamics of the system are determined by a Liouville equation where · · We observe a discontinuity at h = 1. As noted in the text, these transitions are purely nonequilibrium phenomena and they are also present for finite temperatures. Parameters: Therefore, as has been shown in [8,[35][36][37], NESS is a Gaussian state and is fully characterized by the l m l m l m , , of which time evolution is determined by T i r where M r and M i denote real and imaginary parts of the reservoir matrix M, which represents the influence of the environment and depends on the details of the dissipator 4 . The time derivative of the correlation matrix is zero in NESS, hence the steady state correlation matrix C is obtained as a solution of a Lyapunov equation T i The existence of NESS guarantees that equation (16) has a solution, which can be computed efficiently in  N ( ) 3 steps.

Lindblad mesoreservoir
Recently, a model to describe the NESS of the systems, including the degrees of freedom of the finite reservoirs, was introduced [17,38,39]. The idea is not to trace out all the degrees of freedom of the reservoirs but rather to keep some representative parts (mesoreservoirs), which are physically interpreted as the contacts of the system with the reservoirs. The time evolution of the total system (mesoreservoirs and the system of interest) is modeled with the Lindblad master equation such that the mesoreservoirs are thermalized if the couplings between the system and mesoreservoirs are zero. In [17,40] it was shown that in the limit of weak coupling (bright blue) and h = 1 (dark blue). In all cases we observe a logarithmic divergence of the QMI with respect to subsystem size. Other parameters: . 4 The expression for the reservoir matrix M for the mesoreservoir case as well as the derivation of equation (15) are given in appendix A. The derivation and the explicit form of the reservoir matrix M for the Redfield case can be found in [8].
and large mesoreservoir size, the system is thermalized if two mesoreservoirs with the same thermodynamic variables are attached, while the particle current follows the Landauer formula if two reservoirs with different thermodynamic variables are attached. Therefore, one can argue that the mesoreservoir approach gives a meaningful description of a thermal reservoir. In the previous works, particle-conserved Hamiltonians attached to linear dispersion mesoreservoirs have been discussed. Although the linear dispersion mesoreservoir gives the same Liouvillian spectrum (see appendix A), their thermalization properties can be different.
To mimic the * C algebra setup discussed in the previous section, we treat parts of a chain as mesoreservoirs. The total Hamiltonian consists of a one-dimensional XY spin 1/2 chain with a system size = where η k,L and η k,R are diagonal modes of H L and H R , respectively. The operators ν L k, ,1 and ν L k, ,2 can be interpreted as couplings between mesoreservoirs and the traced-out degrees of freedom (super-reservoir).
In subsection 3.1.1 we first study the equilibrium properties of the model. In particular, we will discuss the divergence of magnetic susceptibility 3.1.1. Equilibrium properties of the mesoreservoir approach. In this subsection we first study magnetization and the corresponding susceptibility in the equilibrium state by the Lindblad mesoreservoir approach. We observe a highly fluctuating magnetization profile for < h h c , and a flat magnetization profile except at the boundaries between the system and mesoreservoirs for ) calculated with the mesoreservoir approach agree with the * C results. On the other hand, long-range correlations saturate as a function of | − | l m (for finite Γ and K) to a plateau in the mesoreservoir approach, while they quickly decay exponentially in the * C algebra method. As expected, the agreement between the * C and mesoreservoir correlations is improved for smaller coupling to the super-reservoir Γ (see figure 6(b)) and larger mesoreservoir size K (see figure 7). If the temperature is high enough for low and high temperatures at different values of the magnetic field h. The mesoreservoir results match with the * C algebra for low temperatures, but are notably different for high temperatures. Fortunately, in the regime where we want to observe the QPT, the * C algebra and the Lindblad mesoreservoir approach agree well. Nevertheless, the thermalization at h = 1 is very subtle since the susceptibility diverges logarithmically. To discuss the divergence at h = 1, one should take a limit Γ → 0 and → ∞ K , as can be seen in figure 7, where we show the mesoreservoir size K dependence of the susceptibilities with several coupling strengths Γ at = T 0 L,R . We see that the divergence of susceptibility for is reconstructed with the Lindblad mesoreservoir approach in the limit of Γ → 0 and → ∞ K . The divergence is logarithmic with respect to K and T L,R . Imaginary parts of the correlation matrix vanish linearly with respect to Γ, while specific observables such as heat current vanish for arbitrary Γ.
3.1.2. NESS properties of the mesoreservoir approach. In this subsection, we shall discuss the NESS of the mesoreservoir approach. We numerically observe that the real part of the correlation matrix is an average of equilibrium correlation matrices. Following the discussion in the previous subsection, the divergence of susceptibility is concluded in the mesoreservoir approach. Because of this property, we focus on the imaginary part of the correlation matrix, which represents a genuine nonequilibrium property. Let us first discuss the two extreme cases (Γ ≪ 1 and Γ ≫ 1). For very small Γ the total system (system + reservoirs) is decoupled from super-reservoirs, and the system is expected to follow an average of Fermi distributions element is related to the deviations of the reservoirs' particle occupations from the Fermi distributions [40]). Therefore, one should carefully choose coupling Γ to deal with the imaginary parts of the nonequilibrium case, and in fact, Γ should be in the order of ΔE K / to describe coherent transport in NESS, where ΔE is the width of the mesoreservoirs' energy band. In figures 8(a) and (b) we show comparisons between the NESS correlations of the * C algebra and the mesoreservoir at different magnetic fields (a) and dissipation strengths Γ (b). We observe that the correlations of the * C algebra and the mesoreservoir in NESS quickly differ with the distance from the diagonal elements. Despite the excellent agreement for nearly diagonal elements, we were unable to obtain the jumps in the derivatives of the correlation function. Moreover, we have to choose the right ratio between the coupling strength Γ and the mesoreservoir size K, as can be seen from figures 8(c) and (d), where we show the off-diagonal correlations for different coupling strengths and a fixed mesoreservoir size K = 1600. Since the ratio between Γ and K is important for the nonequilibrium correlations, one sees the strong K dependence on the correlations, contrasting with the fact that real parts of the correlations do not drastically depend on the mesoreservoir size K. For instance, the real parts shown in figure 6 are  quite robust by changing K (thus, we show only K = 400); on the other hand, figure 8 clearly shows that imaginary parts strongly depend on K. Moreover, as discussed in [40], the Γ dependence of is a monotonic decreasing function of K.
In conclusion, similar to the boundary-driven Lindblad model, we also observe signatures of the NQPT in the mesoreservoir approach at = h h c , namely a change in the sensitivity of local observables to the system parameters. Below the critical field = h h c we see large fluctuations, whereas above the critical field there are no fluctuations as we vary the model parameters (see figures 4 and 5). However, these fluctuations are suppressed as Γ goes to zero. On the other hand, we were unable to observe the scaling of the QMI and the far-from-diagonal correlations due to the limitations of the resources, since we should take the limit of large mesoreservoir size and small coupling to the super-reservoir. Nevertheless, the behavior of the off-diagonal correlations shown in figure 8 suggests that the mesoreservoir approach in the limit (from bright to dark) and h = 0.6 (c), h = 0.9 (d). We show the relative values with respect to the appropriate * C algebra result (red dashed lines). In (c), (d) we use K = 400. Other parameters: n = 100 and γ = 0.5 for (a)−(d).
Γ → 0 and → ∞ K resembles the behavior of the * C algebra, i.e. the power law decay of the NESS correlations and the logarithmic divergence of the QMI in all regimes. We also recover the * C algebra divergence of the magnetic susceptibility at zero temperature and h = 1 in equilibrium and nonequilibrium situations. Other genuine nonequilibrium transitions observed with the * C algebra could not be recovered in the mesoreservoir approach.

Modified Redfield master equation
The modified Redfield master equation was studied in [8]   where Γ ω is a Fourier transform of the correlation function Γ μ ν . This choice of correlation functions and coupling operators represents standard reservoirs of harmonic oscillators at two ends with possibly different temperatures. Note that due to the extension of the lower integral bound to minus infinity in the dissipator (19), the frequency cutoff in the spectral functions is irrelevant. In this paper we use the following parameters: In case of equal temperatures and large system sizes 6 → ∞ N ( ), we necessarily recover the results of the * C algebra approach for all values of the anisotropy γ and magnetic field h (see figure 9 (left)). However, in the nonequilibrium setting, where one temperature remains constant and the other goes to zero, the magnetic susceptibility remains finite even for large spin chains and h = 1 (see figure 9 (right)). We also compare the * C algebra and Redfield NESS expectation values of the magnetization for different values of the magnetic field, and observe disagreement around h = 1 (see figure 10 (left)). In figure 10 (right) we show the numerically computed susceptibility. We observe large fluctuations below the critical magnetic field = h h c . These large fluctuations are in agreement with previously observed fluctuations of local observables and the characterization of the long-range magnetic correlation phase with hypersensitivity of local observables on model parameters [8,37]. For a detailed discussion of the NQPT at = h h c in the XY spin 1/2 chain using the master equation approach, see [7,8,37]. Here, we also compare the nonequilibrium susceptibility at h = 1, which, in contrast to the * C result, remains finite even if one temperature of the reservoirs goes to zero (see figures 9 (right) and 10 (right)). Further, we note that the discontinuities of the first and the third derivative of 〈 〉 † f f Im l m at = h h c and h = 1 were not reproduced using the Redfield mesoreservoir approach. Moreover, we find that the imaginary part of the correlation matrix strongly depends on the dissipation strength Γ.

Conclusions and discussion
We studied nonequilibrium quantum phase transitions in the XY spin 1/2 chain analytically using the * C algebra method. First, we showed that the QPT at h = 1 is also present in the nonequilibrium steady state if one temperature of the reservoirs remains absolute zero. In other words, QPT persists even with strong thermal noise coming from one of the reservoirs if the other reservoir is at absolute zero temperature. At the critical point γ = 0 and h = 1, the logarithmic divergence of susceptibility becomes algebraic. Second, we discovered two new transitions which do not exist in an equilibrium state. To be precise, we found discontinuities of the first derivative (at = h h c , γ ≠ 0) and the third derivative (at h = 1, γ ≠ 0) of the imaginary part of the correlation matrix. The former transition appears because a part of the normal modes changes its velocity sign when the magnetic field is smaller than h c , and thus those modes change the information of the reservoirs they carry. On the other hand, the physical interpretation of the latter transition is still unclear. Moreover, at the critical point γ = 0 and h = 1, the jumps in the derivatives disappear. We use these transitions to test the utility of two time generators commonly used in theories of the reduced density operator, namely the Lindblad and Redfield master equations.
• Comparison with the Lindblad mesoreservoir approach. We show that the Lindblad mesoreservoir quantitatively reproduces the QPT in equilibrium. However, we observed that off-diagonal elements which have small expectation values disagree, even in cases of equilibrium. For a nonequilibrium state, we numerically observe that the real part of the correlation is an average of equilibrium values, which agrees with the * C algebra, and does not agree with the modified Redfield equation. For the correlation matrix, we found that the Lindblad mesoreservoir and the * C algebra agree only for correlations near diagonal elements, and are more accurate for small magnetic fields and low temperatures. Despite the good agreement for nearly diagonal elements, we were not able to recover nonequilibrium phase transitions except for the divergence of the susceptibility at h = 1, which seems to be induced by the same mechanism as in an equilibrium case. This may be the effect of finite mesoreservoir size K and coupling strength to environment Γ.
• Comparison with the modified Redfield master equation. The modified Redfield equation by construction exactly describes equilibrium states, which are described by the * C algebra in the thermodynamic limit. Nevertheless, it does not reproduce any nonequilibrium phase transitions observed by the * C algebra. Moreover, we find that the imaginary part of the correlation matrix strongly depends on the dissipation strength Γ.
For the reduced density operator methods (Redfield and Lindblad), the hypersensitivities to the model parameters below h c are reported. In the mesoreservoir case, the fluctuations are suppressed if small dissipator strength Γ and large mesoreservoir size K are taken. Thus, a drastic change in a systemʼs properties at = h h c is a common feature of the * C algebra and the reduced density operator methods, but they are quite different. In the previous works, the transition obtained by a reduced density operator was characterized by the appearance of the correlation resonances [37,41] and different scaling of the QMI in the long-and short-range correlation regimes. On the contrary, it was shown with the * C algebra that the scaling of the correlation elements with the distance from the diagonal remains unchanged as we cross the critical magnetic field h c , i.e. the exponential decay of correlations in the equilibrium and power law decay of correlations in the nonequilibrium case [29]. We also numerically showed that the QMI scales logarithmically with the subsystem size n in all regimes. Therefore, we conclude that the transition obtained by the reduced density operator [7,8] is a consequence of the approach itself.
Since none of the discussed reduced density operator approaches could describe all the transitions obtained by the exact calculations ( * C algebra), the question 'Can any reduced density operator methods thermalize the XY spin 1/2 chain in the complete range of parameters, and at the same time reproduce the nonequilibrium phase transitions obtained by the * C algebra?' remains open. Further, another interesting question arises, namely 'What kind of nonequilibrium quantum phase transitions can we obtain by using different approaches?'. In other words, to what extent do the nonequilibrium quantum phase transitions and the nonequilibrium properties of systems in general depend on the reservoirs (models of open system evolution). transformation Following [29,30], we summarize the setup of the * C algebra approach and the Araki-Jordan-Wigner transformation.
An algebra  is called a * C algebra if it is together with an involution*:   → and finite norm | · |, and satisfies the following properties: , where ᾱ denotes the complex conjugate of α.
•  is complete with respect to a norm | · | < ∞.
The norm-completion of the algebra generated by the Pauli spin matrices forms the * C algebra  S , and the infinite extension of the Hamiltonian (1) defines a * C dynamical system whose dynamics are given by a group of strong continuous *-isomorphism, which is formally given by τ = − A e Ae ( ) t itH itH . In this approach, states of the system ω · ( ) are represented by positive functionals over the * C algebra, and the Hilbert space is introduced as a representation of ( ω , S ). Physically speaking, a * C algebra  represents a set of observables with finite expectation values, a group of strong continuous *-isomorphism τ t gives a time evolution over the * C algebra, and a state ω gives a correspondence between the observables and expectation values , . The equilibrium states with a given τ t and at temperatures T are defined as the states σ · ( ) satisfying the KMS condition: σ τ σ = A B BA ( ( )) ( ).

T i/
As we explained in the main part, a system is initially decomposed into three parts (left semiinfinite, finite system, and right semi-infinite parts). Then, the initial condition is given by where ω ν ν T are the KMS states at temperatures ν T (ν = L, R) of the left and right distinct parts, respectively. The existence of the unique NESS associated to this initial condition ω τ ω = →∞ +