Numerical linked-cluster expansion for the dissipative XYZ model on a triangular lattice

We generalize the numerical linked-cluster expansion (NLCE) method to study the dissipative quantum many-body system. We apply the NLCE to the triangular strip and two-dimensional lattice system. We investigate the dynamics and steady-state properties of the dissipative XYZ model where the coherent dynamics is governed by the anisotropic Heisenberg Hamiltonian while the nonunitary process is induced by the incoherent spin flips. By comparing with the quantum trajectory simulations, the NLCE results show good performance in capturing the dynamics of system with short-range correlations. For strong and long-range correlated system, the larger size clusters in the series should be included. The NLCE study for the magnetic susceptibility also signals the steady-state paramagnetic-ferromagnetic phase transition in the two-dimensional case.


Introduction
The cooperative phenomena in quantum many-body system have attracted much interest in modern physics. It has not only the theoretical significance in the fields of quantum physics, condensed matter and statistical physics, but also the potential applications in the areas of quantum information and quantum optics. Such emerging collective behavior originates from the competitions among the interactions in the system and has been widely studied in the equilibrium case [1]. However, a system is inevitably coupled to the external environment and will be driven away from equilibrium by dissipations. In the presence of dissipations, the evolution of the system is nonunitary and, in the long-time limit, the system will approach to a steady state which is independent of the initial conditions. In contrast to the equilibrium case, the properties of the dissipative system are featured by the evolved steady state rather than the ground state of the Hamiltonian. Recently, a large body of work has be devoted to investigate the properties of the dissipative quantum many-body systems including the dynamics [2][3][4][5][6], transport properties [7][8][9][10][11], steady-state phases [12][13][14][15]16] and phase transitions [17][18][19][20][21], quantum states engineering [22][23][24], as well as the possible applications in quantum sensing [25,26]. Thanks to the experiment progress, dissipative phase transitions have been observed in circuit QED arrays [27] and ensembles of Rydberg atoms [28,29].
Usually the evolution of a dissipative quantum many-body system is governed by the Lindblad master equation under Markovian approximation [30]. The steady state can be obtained by either evolving the master equation for a sufficient long time or diagonalizing exactly the Liouvillian [31]. Because the steady-state solutions in the thermodynamic limit have been proven to be remarkably difficult, some approximations are imposed to the density matrix, such as the single-site and cluster Gutzwiller mean-field factorizations [12,18,16,[32][33][34], to unravel the many-body master equation. Besides, numerical methods based on tensor networks [35][36][37][38], corner space renormalization [39], variational principal [17,40,41], and the neural networks [42][43][44][45] have been proposed.
In this work, we adopt an alternative method, the so-called numerical linked-cluster expansion (NLCE), to study the properties of the dissipative XYZ model in a triangular lattice. On the one hand, the NLCE method enables us to directly access the properties of an infinite lattice system in terms of a sum over contributions of a number of small-size clusters that embedded in the full lattice. It has been used to calculate the equilibrium Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
properties of many-body system [46,47], quench dynamics [48][49][50] and driven-dissipative system [51]. On the other hand the dissipative spin-1/2 XYZ model has attracted a considerable attention. It is a lattice of spins with anisotropic Heisenberg interactions and local incoherent spin relaxations. The competition between the unitary evolution and the incoherent dynamics gives rise to quite interesting phenomena in the dynamics and the steady state. The existence of ferromagnetic steady-state phase which breaks the Z 2 symmetry of the Lindbladian on the two-dimensional square lattice has been predicted [12,18,19,32,33]. Here, we investigate the properties of the dissipative XYZ model on the two-dimensional triangular lattice. The purpose of this work is to explore the potential of NLCE in studying the dynamics of a dissipative many-body system. Although the dissipative manybody systems usually reach the asymptotic steady state in the long-time limit, it is predicted that the limit cycle may exist which breaks the time-translation symmetry of the master equation [16,52]. Additionally the geometry of triangular lattice brings the frustration effect, which refers to competing interactions cannot be satisfied simultaneously, and will enrich the steady-state phases as well as the time-dynamics of the dissipative quantum many-body systems [53,54].
The paper is organized as follows. In section 2, we introduce the dissipative XYZ model and the corresponding master equation. In section 3, we explain the idea of NLCE method. In section 4 we apply the NLCE method to the dissipative XYZ model in the triangular strip and two-dimensional lattice. In section 5 we summarize and present some perspectives.

The model
We consider a spin-1/2 system whose coherent internal dynamics is governed by the following anisotropic Heisenberg Hamiltonian (we will work in units of =  1 ), withŝ a m (α=x, y, z) denoting the Pauli matrices on the m-th site. We assume that the spin-spin coupling is between nearest-neighboring sites. The parameters J x , J y and J z are the coupling strength for the x, y, and z components, respectively. Each spin is subject to an incoherent dissipative process that tends to flip it down along the z-direction at a rate γ.
In the Markovian approximation, the evolution of the dissipative XYZ model can be described by a Lindblad master equation for the density matrixr, [ˆˆ]ˆˆˆ{ˆˆˆ} ( ) å r r r g sr s ss r º =-+ - In the thermodynamic limit, this symmetry may be broken at some region of the parameter space and thus leads a paramagnetic to ferromagnetic steady-state phase transition.
The density matrix at an arbitrary timeˆ( ) r t can be obtained by integrating equation (2). In particular, the steady-state properties of the system can be investigated by looking at the long-time limit (  ¥ t ) of the solution to equation (2) and computing appropriate expectation value of a local observableÂ. An alternative way of finding the steady state sŝ r would be exactly diagonalizing the Liouvillian . The steady state is just the eigenstate of  corresponds to the zero eigenvalue. However, both routes are challenging for a quantum many-body system. On the one hand the dimension of the Hilbert space increases exponentially with the number of sites increasing, on the other hand the incoherent dynamics usually makes the system end up in a mixed state which should be described by the density matrix. Actually, for a N-site system one needs to manipulate a density matrix of dimensions 2 N ×2 N which becomes a computationally intractable task already for a small number of sites (  N 10 ). If one prefers to diagonalizing the Liouvillian, the situation becomes even worse because the dimension of  is the squared ofr. In order to investigate the property of the system in the thermodynamic limit (  ¥ N ), we will adapt the NLCE method in this work. As will be discussed in the next section, the advantage of the NLCE method is that one can access the property in the thermodynamic limit by only manipulating a couple of finite-size systems with small number of sites.

Numerical Linked-cluster expansion
We review the idea of NLCE method in this section. An extensive quantity A of a lattice of size N can always be expressed by a sum over contributions from all clusters that can be embedded in the lattice, where n denotes the size of cluster c n and W A (c n ) is the weight of c n for the property A. The weight is defined according to the inclusion-exclusion principle,  4 . We notice that if cluster c n is formed out of two disconnected clusters c n 1 and c n 2 , i.e. = + c c c n n n . Although the number of clusters in the sum of equation (3) is tremendous, it can be dramatically reduced by grouping together the clusters with common characteristics. If two clusters have the same Hamiltonian, their dynamics should be identical and we need solve the master equation for only one of them. These two clusters are said to be topological equivalent. Group together the topological equivalent clusters, we can rewrite equation (3) as where the inner sum accounts for all topological distinct clusters of size n and l(c n ) is the number of c n appearing in the lattice. The property A per site in the thermodynamic limit can be computed by taking the limit of is the multiplicity of c n presenting the number of embeddings of c n per site. The convergency of equation (6) is guaranteed by the definition of weights. Actually, the averaged A in the thermodynamic limit can be approximated by truncating n=R in the sum which means only clusters at a most size of R are considered. This cutting is reasonable if the typical length scale of correlations in the lattice is less than R.
The primary task in implementing NLCE is to generate all the possible topological distinct clusters for a given size and their multiplicities. An efficient strategy have been explained in [55]. We finish this section by noting that once the topologically distinct clusters and their multiplicities have been generated for a given lattice geometry, the NLCE method can be employed to investigate the properties of systems within in the same geometry with various Hamiltonians and dissipations.

Triangular strip
We start by investigating the dynamics of a strip model consisted of triangles as shown in figure 1(c). The evolution of the density matrix is governed by the master equation expressed in equation (2). We assume that the initial state is a product stateˆ( ) ⨂| | r = +ñ á+ 0 In figure 2, we show the evolution ofˆ( ) s á ñ t z evaluated by NLCE at different orders of R. As a benchmark we perform the quantum trajectory (QT) simulations [56] on a 2×8 strip using open boundary condition, the results are shown by the thick solid lines. In figure 2(a), we show the results for J y /γ=1. One can see that at the initial stage (γ t  0.2), the results for R=5,6,7 converge well because the initial state is uncorrelated. As the time passes, the correlations between distant sites are built which is manifested by the slight deviation of the NLCE results at different orders in the intermediate time (0.2γt0.6). In the long-time limit (γ t  0.6) the curves converge again and reach the steady state. For the chosen parameters the NLCE results at orders of  R 5 already agree with the QT simulation during the whole process of the time-evolution. In this case, within the NLCE method we need only manipulate density matrices with maximal dimension 2 7 ×2 7 . In contrast, in the QT simulation we have to manipulate the state vectors of dimension 2 16 .
However, the power of NLCE method in tracing the dynamics depends on the time-scale over which the correlations are built up. In figure 2(b), we show the evolutions of ( ) s á ñ t z for J y /γ=3 evaluated by the NLCE and QT methods. One see that the converged region of the NLCE results up to R=8 shrinks to γt0.3. Moreover, the converged region extends as R increasing. This tendency implies that one can access the dynamics the by including larger size of clusters in the sum.

Two-dimensional lattice
Now we discuss the application of NLCE in the 2D triangular lattice. The mean-field approximation predicts a paramagnetic to ferromagnetic steady-state phase transition. The steady state in the former phase holds the  2 symmetry while the latter breaks this symmetry by a nonzero magnetic component at the x-y plane. In order to signal the phase transition, we investigate the magnetic linear response to an applied magnetic field in the x-y plane which modifies the Hamiltonian as the following,ˆ(ˆˆ)  , y). We concentrate on the angular averaged magnetic susceptibility, is the induced magnetization, per site, along the direction of the field.
The NLCE study on the steady-state susceptibility is implemented via three steps: (i) compute the induced total steady-state magnetization along α-direction ss of each topological distinct cluster c n ; (ii) compute the per-site α-direction magnetization sŝ s á ñ a in the thermodynamic limit with NLCE; (iii) compute the components c ab of the susceptibility tensor and consequently the averaged av c according to equation (8).
The information of the topological distinct clusters for the two-dimensional triangular lattice have been reported in [47]. We generate these clusters independently and truncated the sum up to a maximal size of R=8.
In figure 3, we show the angular averaged magnetic susceptibility as a function of coupling strength J y /γ at different orders R. For g <  J 0.9 1 y , one can see a good convergency for the present orders, however for J y /γ>1 there is a sudden increase of av c which indicates the PM-FM phase transition. This is consistent with the critical point predicted by the single-site mean-field approximation at J y /γ≈1.017 [12]. The singular behavior of av c in the two-dimensional triangular lattice is similar to that in the two-dimensional square lattice [19,51]. It should be remarked that the NLCE results at different orders for J y /γ>1 has no physical meaning.

Summary
In this paper, we have generalized the NLCE method to the open quantum many-body system. It allows us to access directly the thermodynamic limit and to compute the extensive property of the system. To be concreteness, we applied the NLCE to the dissipative spin XYZ system on the triangular strip and lattice. In the case of a triangular strip we study the time-evolution of the magnetization along z-direction, we find that the NLCE method can perfectly reveal the dynamics with clusters of relative smaller size (R<6) when the correlations in the system is short ranged. However, for a strong correlated system, one need include larger size clusters in the NLCE method to trace the dynamics. In the case of two-dimensional lattice, the NLCE study on the angular averaged magnetic susceptibility converge fairly good far from the critical point. It shows a sudden increasing at J y /γ≈1 indicating the steady-state phase transition from the paramagnetic to ferromagnetic phase.
The present work shows good performance of NLCE method in studying both the dynamics and steady-state property. There are two routes to improve the convergence. The first one is to combine the NLCE with numerical algorithms like corner-space renormalization [39] and Monte Carlo approaches [32] to deal with larger size clusters and include them in the summation. The second route is to utilize some resummation technique, such as Euler transformation and Wynn's algorithm, Brezinski's algorithm [55]. Moreover, other choice of building the clusters in the series expansion may improve the accuracy of the NLCE results [47]. We believe that the NLCE method may be powerful in the general context of the dissipative quantum many-body system. As the future applications of the NLCE method, it will be interesting to check the existence of the limit cycle in the driven-dissipative quantum many-body system which is quite relevant to the concepts of time crystal [57] and quantum synchronization [58]. Moreover, the NLCE studies on the antiferromagnetic steady-state phase in the frustrated dissipative lattice are expected as well.