Nonuniform phases in the geometrically frustrated dissipative XYZ model

The collective behaviors in quantum many-body systems have attracted much attention in modern physics sachdev_book . Usually, frustration leads to exotic properties of many-body system. The frustration refers to the presence of competing interactions that cannot be simultaneously satisfied. A typical example of frustrated system is the spin Ising model on a triangular lattice wannier1950 . Here the frustration stems from the geometry of lattice and thus named as geometric frustration. Although the nearest-neighboring interaction favors the spins being anti-aligned, the three spins on a triangle cannot be antiparallel. Due to the presence of frustration effect, the fluctuations of the ground state are enhanced giving rise to a large degeneracy. A large variety of intriguing quantum phases may emerge in frustrated system, such as spin glass sherrington1975 , spin liquid balents2010 ; zhou2017 , spin ice lee2002 , and supersolid jiang2009 ; zhang2011 . Experimentally, the geometric frustration has been realized in the system of trapped atomic ions kim2010 and optical lattice struck2011 .

However, a system always becomes open due to the unavoidable interactions to its surrounding environment breuer_book . The system-environment coupling leads to the dissipative process of system and drives the system away from equilibrium. In the presence of dissipation, the dynamics of system is described by the master equation. Rather than the ground-state properties in the equilibrium many-body system, we are interested in the long-time limit properties of the open many-body system. In the long-time limit, the system may approach to an asymptotic steady state which corresponds to a fixed point in the phase space. The steady state can break the symmetry of the master equation and constitute novel phases that never appear in the equilibrium counterpart lee2011 ; lee2013 ; jin2013 ; leboite2014 ; weimer2015 ; schiro2016 ; landa2020a . Alternatively, the system may end up in a limit cycle showing a periodic oscillation which intimately relates to the quantum time crystals gong2018 ; iemini2018 ; tucker2018 ; gambetta2019 ; landa2020b and quantum synchronization ludwig2013 ; DT2018 . Recently, the steady-state phase transitions have been experimentally observed in circuit QED lattices fitzpatrick2017 ; collodo2019 and ensemble of Rydberg atoms carr2013 .

Essentially, the peculiar behavior of the steady state originates from the competition between the coherent dynamics governed by the many-body Hamiltonian and the dissipative process. In addition, it is found that the presence of frustration can indeed enrich the steady-state properties. For instance, the geometric frustration in the dissipative spin system leads to the various antiferromagnetism and chaotic dynamics qian2013 ; The steady-state of frustrated driven-dissipative photonic cavities exhibit strongly correlated dark state rota2017a and is analogous to the spin-liquid phase if the two-photon driving and losses are considered rota2019 .

In this work, we study the dissipative spin-1/2 XYZ Heisenberg model on a two-dimensional triangular lattice. By employing the state-of-the art numerical approaches, we explore the steady-state phases of this model. We notice that the same model in two-dimensional square lattice have been widely studied with diverse numerical approaches by focusing only on the steady-state ferromagnetic phase lee2013 ; jin2016 ; rota2017b ; kshetrimayum2017 ; biella2018 ; rota2018 ; nagy2018 ; casteels2018 ; huybrechts2019 . By means of the numerical linked-cluster expansion, the dynamics and phase transition from the paramagnetic to ferromagnetic phases of the dissipative XYZ model on triangular lattice have been investigated qiao2020 . Here, we concentrate on the antiferromagnetic region where the geometric frustration plays crucial role in determining the steady-state phase. We uncover the various antiferromagnetic phases with different configurations of magnetization. Moreover, we also find the evidence of the spin-density-wave (SDW) phase.

The paper is organized as follows: In Sec.II, we explain the frustrated dissipative spin-1/2 XYZ model on triangular lattice and the corresponding master equation that describes the evolution of the system. We also present the possible steady-state phases that may appear in the phase diagram. In Sec. III, as a preliminary exploration, we implement the single-site mean-field (MF) analysis to the steady-state phase. The phases are characterized by the steady-state magnetization which are the stable solutions to the MF master equation. In Sec. IV, we refine the MF results by the cluster mean-field (CMF) approaches in which the short-range correlations are systematically included as the size the of cluster increases. We check the existence of the ferromagnetic, antiferromagnetic and SDW phases in the thermodynamic limit through the Liouvillian gap and the spin-structure factor. We summarize in Sec.V.

The model we consider here is a spin-1/2 XYZ system on a triangular lattice whose Hamiltonian is given by (setting $\hbar$ = 1 hereinafter),

$$H=\sum_{\langle j,k\rangle}(J_{x}\sigma^{x}_{j}\sigma^{x}_{k}+J_{y}\sigma^{y}_{j}\sigma^{y}_{k}+J_{z}\sigma^{z}_{j}\sigma^{z}_{k}),$$

(1)

where $\sigma^{\alpha}_{j}$ $(\alpha=x,y,z)$ are the Pauli matrices on the $j$-th site and the sum of spin-spin interactions runs over the nearest neighboring sites $\langle j,k\rangle$. In addition, we assume that system is coupled to the Markovian environment. Effectively, each spin is subjected to a local dissipative channel that incoherently flips the spin down to $z$-direction. Therefore, the dynamics of the system is governed by the following master equation in the Lindblad form,

$$\frac{d\rho}{dt}={\cal L}[\rho]=-i[H,\rho]+\sum_{j}\mathcal{D}_{j}[\rho],$$

(2)

where ${\cal L}$ is the Liouvillian superoperator and the local Lindbladian superoperator reads

$$\mathcal{D}_{j}[\rho]=\frac{\gamma}{2}\left[2\sigma^{-}_{j}\rho\sigma^{+}_{j}-\{\sigma^{+}_{j}\sigma^{-}_{j},\rho\}\right],$$

(3)

with $\gamma$ denoting the decay rate and the operators $\sigma^{\pm}_{j}=(\sigma^{x}_{j}\pm i\sigma^{y}_{j})/2$ are the raising and lowering operators for the $j$-th spin, respectively. In deriving the Lindblad master equation, the Born-Markovian and the secular approximations are considered. These assumptions can be met in the experimental situations with circuit QED lattice fitzpatrick2017 ; collodo2019 . Recently, other type of master equation beyond Lindblad form have also been developed to study the steady-state properties of open quantum many-body systems lebreuilly2017 ; xu2019 ; reis2020 .

We are interested in the steady state of system which is given by the asymptotic solution to Eq. (2) in the long-time limit $\rho_{\text{ss}}=\lim_{t\rightarrow\infty}{\rho(t)}$ (we use the subscript ‘ss’ to denote steady-state). The steady-state property of the system is characterized by the expectation value of an appropriate operator $\langle O\rangle_{\text{ss}}=\text{tr}(O\rho_{\text{ss}})$.

The master equation in Eq. (2) possesses a $\mathbb{Z}_{2}$ symmetry, which is invariant under a $\pi$-rotation along the $z$-axis $(\sigma^{x}_{i}\to-\sigma^{x}_{i},\sigma^{y}_{i}\rightarrow-\sigma^{y}_{i},\forall i)$. In the thermodynamic limit, as the coupling constants changing, the steady state of the system may undergo a phase transition from the disordered paramagnetic (PM) phase, with vanishing magnetization in the $x$-$y$ plane, to ordered phase, with nonzero $\langle\sigma^{x}\rangle_{\text{ss}}$ and $\langle\sigma^{y}\rangle_{\text{ss}}$, corresponds to the spontaneously broken of the $\mathbb{Z}_{2}$ symmetry. In particular, the ordered phase can be uniform with identical nonzero $\langle\sigma^{x}\rangle_{\text{ss}}$ or $\langle\sigma^{y}\rangle_{\text{ss}}$ through the whole lattice, referred as ferromagnetic (FM) phase. Alternatively, it can be nonuniform with a spatial modulated magnetization of wavelength equal to or greater than the lattice constant, referred as antiferromagnetic phase or SDW phase, respectively. We stress that the anisotropic coupling $J_{x}\neq J_{y}$ is necessary to establish the ordered phases. Because the magnetization along the z-direction is conserved for $J_{x}=J_{y}$. Therefore, there is nothing to counteract the dissipation. The system will eventually approach to the PM phase with all the spin pointing down to the $z$-direction.

The emergence of the ordered phases can be understood as follows. Apart from the dissipation, each spin precesses about an effective magnetic field depending on the direction of its neighboring spin. If the precession is strong enough to counteract the dissipation, the spins will point away from the $z$-direction in the steady state indicating the emergence of ordered phases. Due to the geometric frustration of the triangular lattice, the neighboring spins on the triangles can not align antiparallel in the $x$-$y$ plane in the antiferromagnetic region. This leads to different types of antiferromagnetic orders in steady state as will be discussed in later sections.

We also note that, for $J_{z}=0$, the Hamiltonian in Eq. (1) is reduced to the XY Hamiltonian. The steady-state phases of the dissipative XY model have been reported in Refs. joshi2013 ; wilson2016 ; maghrebi2016 ; marzolino2016 ; parmee2020 . In this paper, we concentrate on the case of $J_{z}/\gamma=1$. The analysis on other values of $J_{z}\neq 0$ works similarly.

We start with the mean-field approximation and will refine the mean-field results by taking the short-range correlations into account in next section. Generally, the mean-field approximation is based on the Gutzwiller factorization by representing the global density matrix as a tensor product of the state of each site, $\rho=\bigotimes_{j}\rho_{j}$. In order to reveal the antiferromagnetic nature we divide the whole lattice into three sublattices $A$, $B$ and $C$ as shown in Fig. 1. Assuming that all the sites belonging to the same sublattice are identical, we may obtain the single-site MF master equation for each sublattice as the followings,

$$\frac{d\rho_{j}}{dt}=-i[H^{\text{MF}}_{j},\rho_{j}]+\frac{\gamma}{2}\left[2\sigma_{j}^{-}\rho_{j}\sigma_{j}^{+}-\{\sigma_{j}^{+}\sigma_{j}^{-},\rho_{j}\}\right],$$

(4)

where $j=A$, $B$ and $C$ denotes the sublattice. The MF Hamiltonian for sublattice $j$ is given by

$$H^{\text{MF}}_{j}=\sum_{\alpha=x,y,z}{\sum_{k\neq j}{J_{\alpha}\langle\sigma^{z}_{k}\rangle\sigma^{\alpha}_{j}}},$$

(5)

where $\langle\sigma^{\alpha}_{k}\rangle$ represent the expectation value of the spin operator for the sublattice $k$. The MF master equations can be reexpressed by the system of Bloch equations as

	$\displaystyle\frac{d\langle\sigma^{x}_{j}\rangle}{dt}$	$\displaystyle=2\sum_{k\neq j}[J_{y}\langle\sigma^{z}_{j}\rangle\langle\sigma^{y}_{k}\rangle-J_{z}\langle\sigma^{y}_{j}\rangle\langle\sigma^{z}_{k}\rangle]-\frac{\gamma}{2}\langle\sigma^{x}_{j}\rangle,$		(6)
	$\displaystyle\frac{d\langle\sigma^{y}_{j}\rangle}{dt}$	$\displaystyle=2\sum_{k\neq j}[J_{z}\langle\sigma^{x}_{j}\rangle\langle\sigma^{z}_{k}\rangle-J_{x}\langle\sigma^{z}_{j}\rangle\langle\sigma^{x}_{k}\rangle]-\frac{\gamma}{2}\langle\sigma^{y}_{j}\rangle,$
	$\displaystyle\frac{d\langle\sigma^{z}_{j}\rangle}{dt}$	$\displaystyle=2\sum_{k\neq j}[J_{x}\langle\sigma^{y}_{j}\rangle\langle\sigma^{x}_{k}\rangle-J_{y}\langle\sigma^{x}_{j}\rangle\langle\sigma^{y}_{k}\rangle]-\gamma(\langle\sigma^{z}_{j}\rangle+1).$

We note that actually there are nine nonlinear differential equations in Eq. (6) since there are three sublattices. For each sublattice, its nearest neighbours consist of the other two sublattices. The steady-state magnetizations of the system are given by the stable fixed points of the nonlinear equations. We numerically find the fixed points by setting Eqs. (6) to be zero and check their stabilities through the eigenvalues of the Jacobian for each fixed point. If the real parts of all the eigenvalues are negative the fixed point is stable, otherwise it is unstable. The presence of two conjugated pure imaginary eigenvalues implies the appearance of limit cycle. In the rest of this paper, we work in units of $\gamma$.

On the triangular lattice, due to the presence of geometric frustration, the antiferromagnetic phase is much richer than the counterpart in the square lattice. We find that the steady-state is possible to enter into two types of antiferromagnetic phase:

(a) Biantiferromagnetic (BAF) phase. The BAF phase corresponds to the case that two sublattices have the same nonzero magnetization in the $x$-$y$ plane but different from the third one, i.e. $\langle\sigma^{\alpha}_{A}\rangle_{\text{ss}}=\langle\sigma^{\alpha}_{B}\rangle_{\text{ss}}\neq\langle\sigma^{\alpha}_{C}\rangle_{\text{ss}}$ ($\alpha=x,y$).

(b) Triantiferromagnetic (TAF) phase. The TAF phase corresponds to the case that the steady-state magnetization of the three sublattices are completely unequal, i.e. $\langle\sigma^{\alpha}_{A}\rangle_{\text{ss}}\neq\langle\sigma^{\alpha}_{B}\rangle_{\text{ss}}\neq\langle\sigma^{\alpha}_{C}\rangle_{\text{ss}}$ ($\alpha=x,y$).

We show the mean-field phase diagram in Fig. 2 which is symmetrized about $J_{x}=J_{y}$. The order parameter is chosen as the steady-state magnetization along the $x$-direction $\langle\sigma^{x}_{j}\rangle^{\text{mf}}_{\text{ss}}$. Basically, there are five steady-state phases in the $J_{x}$-$J_{y}$ plane corresponding to the PM, FM, BAF, TAF, and oscillatory (OSC) phases and two bistable regions with respect to the BAF/PM and BAF/FM phases. In the bistable region, the system evolves to two different phases depending on the initial state.

It is possible to analytically determine the critical points of the PM-FM or PM-TAF phase transitions. For the PM-FM phase transition, the critical points are given

$$J_{x,y}^{c}=J_{z}+\frac{1}{16\mathfrak{z}^{2}(J_{z}-J_{y,x})},$$

(7)

while for the PM-TAF phase transtion, the critical points are given by

$$J_{x,y}^{c}=-2J_{z}-\frac{1}{4\mathfrak{z}^{2}(J_{y,x}+2J_{z})},$$

(8)

where $\mathfrak{z}=6$ is the coordinate number of the triangular lattice.

We deliberately consider a horizontal cut of Fig. 2 at $J_{y}=-1.4$ since this line passes through the typical phases of the present model. The steady-state magnetization is shown in Fig. 3(a). It can be seen that the PM-FM and TAF-FM transitions are continuous and, close to the transition, the growth of the order parameter is described by $\langle\sigma^{x}\rangle_{\text{ss}}^{\text{mf}}\sim(J_{x}-J_{x}^{c})^{1/2}$ where $J_{x}^{c}$ is the transition point. In the PM phase, all spins are uniformly pointing down to the $z$-direction and the steady state is $|\downarrow\downarrow\downarrow...\rangle$. Starting from the PM phase, on the one hand, as the coupling strength $J_{x}$ increasing the system enters to the FM phase. In the FM phase, all the sublattices acquire the same nonzero magnetization along the $x$-$y$ plane. The system is still uniform. On the other hand, as $J_{x}$ decreasing, the steady-state magnetizations of two sublattices deviate from the $|\downarrow\rangle$ state and the magnetization component on the $x$-$y$ plane grow along opposite directions while the third one stays in the $|\downarrow\rangle$ state, indicating that the system enters into the TAF phase.

In an intermediate region $-3.35\leq J_{x}\leq-2.25$, none of the asymptotic steady-state phases are stable and the limit cycle appears. This is supported by the behavior the eigenvalues of the Jacobian for the TAF fixed point. From Fig. 3(g), one can see that, starting from the TAF phase, a pair of complex conjugate eigenvalues cross the imaginary axis (namely two pure imaginary eigenvalues present) at the the transition points $J_{x}=-3.35$ and $-2.25$ signaling the appearance of the limit cycle.

As $J_{x}$ continues to decrease, the TAF phase is replaced by the BAF phase at $J_{x}\leq-3.49$ and the magnetization of all the three sublattices have nonzero component in the equatorial plane. The continuity of the PM-FM and PM-TAF transitions can be seen from the value of $\langle\sigma^{x}\rangle^{\text{mf}}_{\text{ss}}$ in the vicinity of the critical points. Additionally, the transition between the nonuniform BAF and TAF phases is discontinuous characterized by sudden change of the order parameter.

In order to give an intuitive description, we show the time-evolution of the magnetization on the equatorial plane of Bloch sphere for each phase in Fig. 3(b)-(f). We initialize the sublattice by the $120^{\circ}$ state on the equatorial plane of the Bloch sphere, i.e. $|\psi_{A}(0)\rangle=(|\uparrow\rangle+|\downarrow\rangle)/\sqrt{2}$, $|\psi_{B}(0)\rangle=(|\uparrow\rangle+e^{i2\pi/3}|\downarrow\rangle)/\sqrt{2}$ and $|\psi_{C}(0)\rangle=(|\uparrow\rangle+e^{-i2\pi/3}|\downarrow\rangle)/\sqrt{2}$. It is found that the system always falls onto the stable fixed points. For the uniform phases, there is only one fixed point. More specifically, for the PM phase the fixed point locates on the origin while for the FM phase the fixed point deviates from the origin, implying the preservation and breakdown of the $\mathbb{Z}_{2}$-symmetry. For the nonuniform phases, there are two and three fixed points for the BAF and TAF phases, respectively. Interestingly, there exists an OSC phase in which the long-time limit trajectory is a limit cycle breaking the time-translation symmetry. The presence of the OSC phase in the square lattice has been ruled out in the literature owen2018 . Here, the presence of oscillatory phase in frustrated lattice is probably an artifact of the sublattice factorization of the full lattice. As will be seen in next section, the OSC region is replaced by the time-translation invariant phases in the thermodynamic limit.

In this section we will investigate the steady-state phases beyond mean-field approximation. It has been proven that the short-range correlations play crucial role in determining the structure of the steady-state phase diagram. Here, we systematically include the short-range correlations in the analysis by exploiting the CMF approach. The idea of CMF approach is to decouple the Hamiltonian of the full lattice into Hamiltonians of identical clusters consists of contiguous sites. The interactions inside the cluster are faithfully described while the interactions between sites belonging neighbouring clusters are effectively treated by mean-field terms. More details about the CMF approach can be found in Ref. jin2016 .

For the present model, we perform the CMF studies with a series of triangular clusters up to size $L=15$ as shown in Fig. 1. For $L=3,6$ and $10$ we simulate the evolution of the CMF master equation with a fourth-order Runge-Kutta algorithm while for $L=15$ we resort to the quantum trajectory. In this work, the data obtained by quantum trajectory are averaging over $500$ trajectories.

In summary, we have investigated the steady-state phase diagram of dissipative spin-1/2 XYZ model on a triangular lattice. We employed the mean-field type approximations based on the Gutzwiller factorization of the global density matrix to decouple the full master equation. In the single-site MF analysis, we numerically solved self-consistent Bloch equations and check the stabilities for the fixed points. Due to the geometric frustrations in the lattice, we found that the MF steady-state phase diagram is rather rich including the PM, FM, BAF, TAF and OSC phases.

In order to confirm the existence of the various steady-state phases in the thermodynamic limit, we gradually included the short-range correlations in the analysis by enlarging the sizes of clusters. We demonstrated that FM phase shrinks into a closed region which is supported by the closing of Liouvillian gap in the $\mathbb{Z}_{2}$ symmetry-broken phase. We also showed that the antiferromagnetic region in the mean-field approximation changed dramatically: first the OSC phase is erased and replaced by the TAF phases; second a new region of BAF emerges in an intermediate region inside the TAF phase. Disappearance of the limit cycles, when the short-range correlations are considered in the analysis, indicates that the OSC phase is an artifact of the MF treatment. Moreover, we found the evidence of the SDW phase, which is missed by the MF approximation, by performing the linear stability analysis. The feature of SDW phase was revealed by the spin-structure factor of a finite-size cluster. To further characterize the steady-state phase transitions, the quantum fisher information is proved to be a potential indicator braun2018 ; marzolino2017 .

Finally, we note that identifying the steady-state phases of open quantum many-body system is rather a demanding task. Due to the exponential growth of the many-body Hilbert space, we stopped our CMF computation at the cluster with $L=15$. To exactly locate the boundary of phases and extract the critical exponents, further analysis on larger size clusters are required. Along this line, the combination of the CMF approximation with other algorithms, such as corner-space renormalization finazzi2015 , matrix-product-operator approach mascarenhas2015 , and neural networks yoshioka2019 ; nagy2019 ; hartmann2019 ; vicentini2019 should be helpful.

This work is supported by National Natural Science Foundation of China under Grant No. 11975064 and No. 11775040, and the Fundamental Research Funds for the Central Universities No. DUT19LK13.