Finite-temperature symmetric tensor network for spin-1/2 Heisenberg antiferromagnets on the square lattice

Within the tensor network framework, the (positive) thermal density operator can be approximated by a double layer of infinite Projected Entangled Pair Operator (iPEPO) coupled via ancilla degrees of freedom. To investigate the thermal properties of the spin-1/2 Heisenberg model on the square lattice, we introduce a family of fully spin-$SU(2)$ and lattice-$C_{4v}$ symmetric on-site tensors (of bond dimensions $D=4$ or $D=7$) and a plaquette-based Trotter-Suzuki decomposition of the imaginary-time evolution operator. A variational optimization is performed on the plaquettes, using a full (for $D=4$) or simple (for $D=7$) environment obtained from the single-site Corner Transfer Matrix Renormalization Group fixed point. The method is benchmarked by a comparison to quantum Monte Carlo in the thermodynamic limit. Although the iPEPO spin correlation length starts to deviate from the exact exponential growth for inverse-temperature $\beta \gtrsim 2$, the behavior of various observables turns out to be quite accurate once plotted w.r.t the inverse correlation length. We also find that a direct $T=0$ variational energy optimization provides results in full agreement with the $\beta\rightarrow\infty$ limit of finite-temperature data, hence validating the imaginary-time evolution procedure. Extension of the method to frustrated models is described and preliminary results are shown.


Introduction
The spin-1/2 Heisenberg antiferromagnet on the square lattice orders magnetically at zero temperature but, compared to the classical Néel state, quantum fluctuations induce a reduction of the order parameter. The value of the staggered magnetization is about 60% of its classical counterpart (1/2) [1][2][3]. At any finite temperature, Mermin-Wagner theorem implies the restoration of the full continuous spin-SU(2) symmetry but the magnetic correlation length diverges exponentially fast when approaching zero temperature [4]. These features have been well established by quantum Monte Carlo (QMC) simulations [1,3,5,6] which are free of any sign problem in the absence of magnetic frustration.
In the last decade, tensor network methods have brought considerable new insights in the understanding of models of interacting quantum spins. At zero temperature, two-dimensionnal Projected Entangled Pair States (PEPS) have provided many valuable ansatze for groundstates satisfying the area-law of the entanglement entropy [7][8][9][10][11]. Although small deviations from the area-law are expected [12,13], simple PEPS with small bond dimensions offer good approximate realizations of the symmetry-broken Néel state [14,15]. However, describing SU(2)-symmetric groundstates with critical or even very long-range correlations remains an open challenge [16].
The implementation of finite-temperature tensor network algorithms in two-dimensions have been steadily developed and improved over the last ten years [17][18][19][20], and benchmarked on simple models as the quantum Ising model [18,19,21,22], the quantum compass model [23] or the Shastry-Sutherland Heisenberg model [24,25]. These methods are based on a smart parametrization of the thermal density operator (TDO) in terms of a double-layer Projected Entangled Pair Operator (PEPO), which naturally 1 arXiv:2010.07828v1 [cond-mat.str-el] 15 Oct 2020 guarantees its positivity. This formalism also allows to construct the TDO via an imaginary-time evolution [26][27][28][29]. The method can be used directly in the thermodynamic limit with infinite-PEPO (iPEPO). Its efficiency is remarkably good, provided the phase is gapped. Here, our aim is to test this method in the more complex case of a gapless phase, with a diverging correlation length at low temperature. In addition, we would like to implement the lattice (translation and C 4v point groups) and spin-SU(2) symmetries within the ansatz itself. 2 Symmetric iPEPO method

Thermal density operator framework
The method is based on a simple tensor network approximation of the (unnormalized) thermal density operator ρ(β) = exp (−βH) at temperature T = 1/β, as depicted in Figs. 1. The building block of this construction is the rank-6 on-site tensor A of Fig. 1(a) containing 4 virtual bonds of dimension Table 1: Numbers of elementary tensors in each of the tensor classes defined by their virtual spaces V. The number of elementary C 4v -symmetric tensors are shown in the third column. When the site C 4v point group symmetry is broken to C s , the N D tensors are "split" into N D tensors (fourth column). For D = 7, assuming "color" spin-1 exchange symmetry reduces the number of elementary tensors, as shown by the numbers in parenthesis. (31) 115 (73) D to be connected to the neighboring sites, a physical bond (in blue) carrying the d = 2 spin-1/2 degrees of freedom and an ancilla bond (in red) of same dimension. The square lattice translation and point group symmetry can be enforced by assuming identical C 4v -symmetric tensors on every site. In addition, since spin SU(2)-symmetry is preserved at any finite temperature, we also enforce the SU(2)-symmetry at the level of the on-site tensor (see later for details). An infinite Projected Entangled Pair Operator (iPEPO) is obtained by contracting the tensor network over all virtual degrees of freedom 1 . The thermal density operator can then be approximated by a double layer of such iPEPOs, after contracting over all ancillas as shown in Fig. 1(b). It has been pointed out that such a construction naturally guarantees the positivity of the (approximate) TDO. Note that this approximative representation of the TDO becomes more and more accurate for increasingly larger (virtual) bond dimension D. With the iPEPO approximation of the TDO at hand, any thermal average can, in principle, be computed by applying a local operator (or a string of operators) on the (e.g. top) physical layer and tracing over all top and bottom physical degrees of freedom.
Symmetric PEPO ansatz -Here, by enforcing the lattice and spin symmetries of the problem, we simplify the ansatz further since the on-site tensor can then be written as a linear superposition of a small number N D of (fixed) symmetric tensors T α , where c α (β) are real temperature-dependent coefficients. The temperature dependence of the TDO ρ β is therefore only encoded in the temperature dependence of the latter coefficients. The construction of symmetric PEPO is a generalization of the standard method introduced to derive symmetric Projected Entangled Pair State (PEPS) tensors [30]. In fact, as depicted in Fig. 1(c), the family of tensors T α can be obtained by contracting two families of PEPS tensors over their "physical" space 0 ⊕ 1. The family of rank-5 A tensors is constructed assuming spin-SU(2) and C 4v symmetry w.r.t. exchange of the four virtual bonds (thick segments). There are two rank-three tensors corresponding to each fusion outcomes -spin-0 or spin-1 -of the physical and ancilla degrees of freedom, and therefore two independent subclasses of A tensors which can be contracted with each of them.
The families of A tensors are defined by the vector space V of the virtual degrees of freedom written as a direct sum of SU(2) irreducible representations (irreps). At infinite temperature, i.e. β = 0, the unnormalized TDO reduces to the identity and the virtual space reduces to a single singlet state, V 0 = 0, D = 1. At finite β we increase the entanglement between sites by introducing more . The number of C 4v SU(2) symmetric tensors T α in each class is shown in the third column of Table 1.

Plaquette Trotter-Suzuki decomposition
With the framework established, one needs a reliable method to obtain the β-dependence of the coefficients c α (β). We shall use here a Trotter-Suzuki decomposition. The usual initial starting point consists in splitting the unnormalized TDO into a large number N τ of imaginary-time evolutionary steps, ρ(β) = [exp (−τ H)] Nτ , where τ = β/N τ is a small imaginary-time step, and start from the infinite-temperature i.e. β = 0 limit. In our framework, assuming that, at a given temperature T = 1/β, ρ β is in the (approximate) iPEPO form, the problem reduces to evaluate the action of the"infinitesimal" propagator G(τ ) = exp(−τ H), i.e. ρ(β + τ ) = G(τ ) ρ(β), and to approximate the result as a new iPEPO defined by an updated tensor A(β + τ ) corresponding to the set of coefficients {c α (β + τ )}. To realize this, in practice, one can rewrite the Hamiltonian as Trotter errors.

Variational optimization of the plaquette fidelity 2.3.1 First step: a plaquettes
Let us first concentrate on the first product term a G a (τ ) of G(τ ). Since all G a plaquette propagators commute with each other, it is sufficient to consider the action of one of them to obtain the new tensors on a single plaquette, and then update the tensors on the whole lattice. Following the notation of Fig. 3(a,b), we can rewrite ρ(β) = LL T and G(τ ) = G 2 . Up to a cyclic permutation of the operators (which does not change thermal expectation values) we obtain ρ(β + τ ) = (LG)(LG) T which we have to approximate as ρ(β + τ ) L L T , where L is a single layer iPEPO which differs from L only on the four sites of a single plaquette, by a new A tensor. Since the action of a G a (τ ) breaks the square lattice C 4v symmetry into C s (only the reflections w.r.t the a plaquette diagonals remain), the A tensor belongs to an enlarged class, where the new class dimensions N D are listed in the fourth column of Table 1. Outside of the active 2 × 2 plaquette, all fidelities involve the same uniform tensor network of onsite C 4v -symmetric double-layer AA T tensor contracted over both physical and ancilla degrees of freedom. We have used a single-site (symmetric) Corner Transfer Matrix Renormalization Group (CTMRG) [30,31] to contract the network around the active plaquette, resulting into a converged (socalled "fixed-point") SU(2)-symmetric environment of adjustable bond dimension χ, as shown e.g. in Fig. 4(a).
Minimizing the distance G|L − |L between G|L and |L is equivalent to maximize the (normalized) fidelity where the overlaps ω 0 , ω 1 and ω 2 are defined as and graphically represented in Fig. 4. We have used a Conjuguate Gradient (CG) optimization routine to maximize O w.r.t. the set of coefficients {c α (β + τ )} defining the A tensor located on the four sites of the 2 × 2 plaquette of |L . Figure 4: (a) Graphical representation of the overlaps ω 0 , ω 1 , and ω 2 , obtained by using the CTMRG fixed-point environment, and by contracting the operators Ω 0 (b), Ω 1 (c) and Ω 2 (d), respectively, over the ancilla degrees of freedom and drawn in (a) using the same generic green box (e).

2nd step: b plaquettes
Once the optimum A tensor is obtained on the 2 × 2 plaquette, it is updated on all the other lattice sites. We then have to consider the new fidelity O defined on a 2 × 2 b plaquette by replacing A → A and A → A , in the bottom and top layers respectively, and optimize it w.r.t. the new C s -symmetric tensor A . Applying G on a b plaquette, instead of an a plaquette, amounts in fact to rotate the four (bottom) A tensors by 180 o . We know that, in the limit τ → 0, the optimum A tensor should become exactly C 4v -symmetric, and the small deviation from C 4v -symmetry is due to the non-commutativity of G a and G b on neighboring plaquettes and, hence is part of the Trotter error. We can partially correct it by symmetrizing the A tensor and update it on the whole lattice for the next Trotter step.

Full environment versus simplified environment
At this stage, it is useful to specify how the overlaps ω i are computed. As mentioned above, we use a single-site CTMRG algorithm to contract the network around the active plaquette. This procedure is parametrized by the number χ of (virtual) states kept at each step of truncation of the corner transfer matrix (see Ref. [30,31] for details). In the course of the imaginary-time evolution (TE), we have used a fixed χ = χ TE , typically set to D 2 or 2D 2 , to compute with CTMRG the new converged C and T environment tensors (see Fig. 4(a)) at each Trotter step. For D = 4, we have used a full environment (FE) scheme, meaning that the environment used to construct the ω i overlaps retain all the χ 2 and χ 2 D 2 components of the fixed-point C and T tensors. In contrast, for D = 7, we used a simplified environment (SE) where we used only the largest weight of the fixed-point χ × χ corner matrix (and the corresponding D 2 components of the T tensor) to construct the environment. Note however that, in both cases, i) the CTMRG involves all χ degrees of freedom on the bonds (in this sense the SE scheme is more elaborate than the simple update scheme [32] for which no CTMRG is performed) and ii) the same type of "brute force" variational optimization w.r.t. the components of the A (or A ) tensor is performed via a CG algorithm, necessitating repeated computations (a few thousand times at each Trotter step) of the overlaps ω i . Also, once the A and A tensors have been optimized, all observables (like the energy density) are computed using the full environment.

Heisenberg model and comparison to QMC
We now focus on the spin-1/2 Heisenberg model on the square lattice, involving only an antiferromagnetic coupling J 1 (set to 1) between nearest-neighbor sites i, j , enabling a direct comparison to quantum Monte Carlo (QMC) results. In that case, the simplest a-b Trotter-Suzuki decomposition of Fig. 2(a) is possible, since all the Hamiltonian terms are contained within the a and b plaquettes. The QMC computations, which are statistically exact, are performed at finite temperature on finite albeit large samples of linear size L, thanks to the efficient loop algorithm [33]. When the finite sample of size L is several times larger than the correlation length ξ (one usually considers L > 6ξ), the QMC results can be considered to be in the thermodynamic limit, thus allowing to benchmark the iPEPO results. QMC simulations are performed on systems of size up to L = 256 (with N = L 2 spins), leading to the absence of finite-size effects down to inverse temperature β = 4 (where ξ 40). The correlation length ξ is estimated from the second moment estimator [34,35] using a loop improved estimator [33,35].

Thermal energy density
As a first benchmark we computed the temperature dependence of the mean energy (per site), e = H T /N where the thermal average is defined as usual as · · · T = Tr(ρ β · · · )/Tr(ρ β ), using the same χ = χ TE as for the overlaps ω i . The results are shown in Fig. 5(a) for the two ansätze corresponding to D = 4 and D = 7, with small time steps for which Trotter errors are almost negligible. At high temperatures, say for β < 1.5, the iPEPO data follow very closely the QMC data 2 . Note that, already for β > 0.5 strong deviations occur w.r.t second order high-temperature expansion, , which is picked-up relatively well by our method. This is clear from Fig. 5(b) where the first-order, linear in β, contribution has been subtracted off to zoom in on the small energy difference: above β ∼ 0.5 the iPEPO data deviate strongly from the second-order prediction while staying quite close to the exact QMC data. We believe that the small remaining deviations still observed for D = 7 in the range 0.5 ≤ β ≤ 1.5 are due to the approximate SE procedure used in that case to compute the fidelity. When lowering the temperature further, the energy density starts to saturate above the exact ground-state energy e QMC 0 −0.6694 (see Ref. [3]), around −0.636 and −0.661 for D = 4 and D = 7, respectively. Note that the ansatz with 0 ⊕ 1 2 ⊕ 1 (D = 6) is not providing a significant improvement compared to D = 4, so we shall not consider it further.

Spin correlation function and correlation length
At any finite temperature, the spin-spin correlations are short-range (from Mermin-Wagner theorem) but it is known that the correlation length ξ associated to the long-distance exponential decay rises exponentially fast at (not so) large inverse-temperature β. For example, at β = 5 (respectively β = 6), ξ is already beyond 120 (respectively 400) lattice spacings [4][5][6]. To investigate how good this feature is captured by our iPEPO approach we have computed the spin-spin correlation function C(|i − j|) = S i · S j versus distance r ij = |i − j|, and extracted the correlation length ξ from an exponential fit exp (−r ij /ξ) of the long distance behavior 3 . Although the tensor optimization is always realized  As shown in Fig. 6(c), showing the scaling of the correlation length ξ with the environment dimension χ, larger values of χ are needed for convergence at larger β and at larger D as well. At high temperature, the iPEPO correlation length reported in Fig. 6(d) increases linearly with β, tracking very well the QMC results. However, above β ∼ 2, the QMC correlation length shoots up while the D = 4 iPEPO correlation length starts to saturate around 7 lattice spacings. Even for larger D = 7, deviations from QMC still occur quickly at intermediate β. For example, while the (exact) QMC correlation length at β = 3 (β = 4) is around 12 (39) lattice spacings, we get ξ 4.2 (ξ 5.6) and ξ 5.5 (ξ 9.2) for D = 4 and D = 7, respectively.
The iPEPO spin correlations at short and intermediate distances have been directly compared to QMC results, obtained on system sizes for which finite size effects are negligible. As shown in Fig. 6(a), a good agreement with the D = 4 results is obtained up to β = 1.5, while significant deviations occur at β = 2. Even for D = 7, strong deviations occur above β ∼ 2, as seen in Fig. 6(b) comparing results at β = 3 and β = 4. Interestingly, we observe that QMC results at β = 3 agree with iPEPO results at β = 4. This suggests that one may define an effective inverse-temperature β eff (β, D) < β such that local observables computed at β approximately match their exact values at β eff . Of course, this would imply an equivalence between the two TDO (not only between specific observables) which is difficult to prove.

Finite correlation length scaling
Although the iPEPO correlation length deviates strongly from the exact behavior at intermediate and low temperatures, qualitative agreement between iPEPO and QMC still survives for some observables providing one considers their behaviors w.r.t. the inverse correlation length rather than w.r.t. the inverse-temperature. Such inverse correlation length scaling has been introduced as a powerful tool for zero-temperature symmetry-broken groundstates [36,37] and used recently to extract quantitatively the staggered magnetization curve of the frustrated quantum Heisenberg model [15]. Also, it has been applied to the case of second-order finite-temperature phase transitions associated to the spontaneous breaking of a discrete symmetry, as in the 2D quantum Ising model [38].
We investigate here whether the same concept still holds in the more challenging case of a zerotemperature phase transition associated to a continuous-symmetry breaking. Fig. 7(a) shows the behavior of the mean energy versus ξ-inverse where different temperature ranges have been considered for iPEPO and QMC, β ∈ [1,8] and β ∈ [1, 4], respectively. The QMC energy varies almost linearly with ξ −1 in rough agreement with the iPEPO data. Note that, although deviations do occur, such a qualitative agreement is nevertheless rather unexpected considering the very important mismatch of the correlation lengths at intermediate temperature.
We have also computed an estimator of the (zero-temperature) staggered magnetization m q = | 1 N i exp(iq · r i ) S i 0 |, with q = (π, π), in terms of the spin structure factor S q = 1 N ij exp(iq · r ij ) S i · S j = r = (−1) r S 0 · S r . When β → ∞, S q should diverge as ξ 2 i.e. as the area of (antiferromagnetically) correlated spins. More precisely, if we crudely assume that spins are correlated, i.e. S 0 · S r (−1) r (m q ) 2 , within a disk of radius ξ and that S 0 · S r 0 beyond, we then obtain when ξ 1. The quantity m q = S q /π 1/2 ξ on the right-hand side of (8) gives then a finite-temperature estimator of the (zero-temperature) staggered magnetization m q . Note that, because of the abrupt separation assumed between correlated and uncorrelated spins, the r.h.s is approximate, even in the β → ∞ limit, and the exact relation between the two quantities should involve a numerical prefactor (close to 1). We have plotted m q versus 1/ξ in Fig. 7(b) and compared it to the QMC estimate. Although the agreement is still rough, our data extrapolated to ξ −1 = 0 clearly indicates that the zero-temperature staggered magnetization is indeed finite.

Zero-temperature limit
So far the optimization of the finite-T site tensor A has been performed via imaginary-time propagation. An alternative method would involve the direct variational optimization of the free energy (per site) f = e − T s, where s is the thermodynamic entropy per site. Computing s is a notoriously hard problem for tensor networks but, at T = 0, it becomes feasible to minimize the thermal energy e = 1 N Tr(Hρ(A)) which only involves the local tensor A to optimize upon. We have accomplished this task using the same global optimization scheme as used for optimizing iPEPS [16,39]. Our T = 0 result for D = 4 reveals a quite large spin correlation length, around 16 lattice spacings, as shown in Fig. 6(c). For D = 7, our CTMRG failed to produce a SU(2)-invariant boundary, probably because the correlation length becomes very large in the course of the optimization inducing spontaneous symmetry breaking. We have also included the T = 0 data point obtained for D = 4 to Figs. 7(a) showing that the finite-T data seem to approach, for larger and larger β, the T = 0 "end-point". This provides a check that our imaginary-time propagation method remains reliable at low temperature in producing the optimal TDO within its variational manifold.
Note that, in Figs. 7(a,b), we have not shown any QMC data beyond β = 4, value of β at which finite size effects start to come into play on the largest 256 × 256 cluster. In the regime β 1, one expects the energy to behave as e = e 0 + Cβ 3 [3][4][5] and then to approach the T = 0 limit e 0 as −1/ ln 3 (ξ −1 ) i.e. with a vertical slope in Fig. 7(a).

Extention to frustrated models
Lastly, we introduce frustrating interactions. The Hamiltonian takes the form : where the antiferromagnetic J 2 interaction couples all next-nearest-neighbor sites. Due to the frustration, quantum Monte Carlo suffers from the sign problem and cannot provide precise results here. Recent Variational Monte Carlo [40], DMRG [41,42], finite PEPS [43] and iPEPS [15] calculations show that the (T = 0) Néel phase survives up to J 2 /J 1 0.46 (2), where a second-order phase transition to a (possibly critical) spin liquid takes place [16,[41][42][43]. We choose here J 2 = 0.5 (J 1 is set to 1 as before), inside the spin liquid phase.
Since interactions occur on all plaquette diagonals, we now have to use the a-b-c-d TS decomposition of Fig. 2(b) 4 , i.e. write the infinitesimal imaginary-time propagator as G(2τ where the J 1 couplings have been equally split between adjacent plaquettes. Then, in addition to the previous two first steps associated to the a and b plaquettes, we now have to complete two more steps. After the two first steps, the optimized A tensor is approximately C 2v -symmetric (instead of C 4vsymmetric, because of the extra diagonal interactions) and we (partially) correct the deviation due to the Trotter error by symmetrizing it. For the c plaquette, the (C 2v -symmetrized) A tensor is rotated by 90-degrees and the optimization of the fidelity leads to a new A tensor, which is rotated by 180degrees and used for the d plaquette. We observed that the last optimized tensor after these 4 steps is almost C 4v -symmetric, as expected, and can be used (after exact symmetrization) in the single site CTMRG to compute the new environment.
The temperature-dependance of the thermal energy versus the inverse-temperature β is shown in Fig. 8 using a Trotter step τ = 1/16. Furthermore, from the series expansion in β up to 9th-order [44], we constructed several Padé approximants which show convergence up to β 2 (we display the specific Pade [4,5] approximant in Fig. 8). Even at the smallest β, for which the energy is linear in β, we see deviations from the expected analytical predictions. This suggests that our iPEPO ansatz is not optimal to approximate the TDO of the frustrated model. We believe, either the dimension D = 7 of the virtual space of the iPEPO is still not large enough or the analytical iPEPO form of the TDO itself is not appropriate. In any case, we observe that, at large β, the thermal energy seems to slowly decrease towards the estimated groundstate energy [16,41] shown on the plot, rather than to exponentially saturate to a minimum value. This may be a signature of the critical nature of the spin liquid phase [16,[41][42][43].

Conclusion and perspectives
Within the last ten years, tensor network methods have become the state-of-the-art technique to deal with two-dimensional correlated lattice models, and frustrated quantum spins in particular, for which QMC is inapplicable due to the sign problem. Here, we have more specifically examined the efficiency and the accuracy of tensor network techniques to compute finite-temperature properties in quantum spin systems with diverging correlation length at low temperature, a particularly challenging situation. In addition, our goal was to explicitly preserve the symmetries of the problem, namely the lattice symmetry as well as the continuous spin-rotation (SU(2)) symmetry. Thirdly, we wished to design a set up capable to deal with longer-range (frustrating) interactions. Our symmetric iPEPO ansatz of the TDO, associated to a plaquette TS discretized time-evolution fulfils such requirements. This framework has first been tested on the square lattice spin-1/2 Heisenberg model and confronted to large-scale QMC results (in the thermodynamic limit). First, we found that our ansatz, despite its relatively small bond dimension, is capable to generate thermal states with large correlation length, of the order of 20 lattice spacings or even more. Unfortunately, the rapid exponential growth of the Heisenberg model correlation length at low temperature is beyond the current ability of the method. Nevertheless, we showed that some observables (like the thermal energy or the staggered magnetization estimator) follow an approximate inverse-correlation length scaling, which could be of practical use for future investigations.
Our framework based on plaquette Trotterisation, allowed to also include antiferomagnetic coupling between next-nearest neighbor sites, and to investigate the long-time debated J 1 − J 2 quantum spin model on the square lattice. Our preliminary results show a "proof of principle" that the technique can be used in such a case. However, we found that our iPEPO does not accurately capture the correct temperature behavior. We believe a more entangled ansatz with increased bond dimension becomes necessary in the presence of magnetic frustration. This challenging issue is left for future studies.