---
abstract: 'Classical frustrated spin systems give rise to many fascinating many-body phenomena, but standard computational techniques such as Monte-Carlo sampling or tensor networks have great difficulties in simulating a large variety of them. We propose a framework to map generic frustrated spin systems onto frustration-free models on a correlated phase space, resulting both in an improved analytical understanding of the ground state ensembles via local rules, and a natural construction of a tensor network on which standard algorithms can be used efficiently. The main technical ingredient consists in a linear program identifying large scale degrees of freedom for which the frustration can be relaxed. We illustrate the power of the method by determining the ground-state local rule and computing the residual entropy of a frustrated Ising spin system on the kagome lattice with next-next-nearest neighbour interactions.'
author:
- Bram Vanhecke
- Jeanne Colbois
- Laurens Vanderstraeten
- Frédéric Mila
- Frank Verstraete
bibliography:
- 'bibliography.bib'
title: Relaxing Frustration in Classical Spin Systems
---

***Introduction.*** One of the most beautiful manifestations of emergent behaviour in statistical physics can be found in the arena of frustrated spin systems [@Lacroix2011]. Frustration in a classical spin system occurs whenever it is impossible to find a spin configuration which minimizes the energy of all terms in the Hamiltonian simultaneously. This phenomenon can lead to macroscopic ground state degeneracies, giving rise to interesting zero-temperature physics such as effective realizations of gauge theories [@Henley2010].

Early exact results in this context were obtained for frustrated Ising models on all planar two-dimensional lattices with nearest-neighbour interactions using a mapping to free fermions [@Kasteleyn1963; @Schultz1964], and for more general systems such as planar spin ice [@Lieb1967] using Bethe ansatz techniques [@Baxter1982]. It has, however, proven difficult to treat frustration in generic (i.e., non-integrable) models: to reach the low-energy phase space and sample it efficiently, Monte Carlo methods require ad-hoc non-local cluster updates to fight both critical slowing down [@Swendsen1987; @Wolff1989] and frustration [@Wang2012; @Rakala2017]. In addition, calculating the free energy requires the use of thermodynamic integration, making zero-temperature residual entropies hard to determine accurately.

Tensor networks [@Verstraete2008] provide a new computational approach for tackling classical lattice models with strong correlations down to zero temperature, as was recently demonstrated by the determination of the residual entropy of ice and dimer configurations in three-dimensional lattices with unprecedented precision [@Vanderstraeten2018]. This was achieved by employing matrix product state (MPS) and projected entangled-pair state (PEPS) algorithms — originally devised for finding ground states of strongly correlated quantum many-body systems — to determine the leading eigenvectors of row-to-row or plane-to-plane transfer matrices. In these applications, the transfer matrix represents an explicit encoding of the local rule for ice or dimer configurations [@Haegeman2016].

However, if a local rule is not known, the application of tensor networks to low-temperature frustrated spin systems is less straightforward. Indeed, in the standard way for representing a partition function, the spins are represented by delta tensors, whereas interactions can be encoded in simple matrices $t$ on the bonds between the spins. The partition function for the triangular-lattice Ising model, for example, is represented as $$\label{eq:Z}
\mathcal{Z} = \diagram{1}$$ where the matrix $t$ contains the Boltzmann weights for the nearest-neighbour interactions $$t_{ij} = \e^{ -\beta \mathcal{H}(s_i,s_j) } = \begin{pmatrix} \e^{-\beta J} & \e^{+\beta J}  \\ \e^{+\beta J} & \e^{-\beta J} \end{pmatrix}.$$ In the zero-temperature limit, we can regularize the partition function as follows: $$\mathcal{Z}_0=\lim\limits_{\beta\to\infty}e^{\beta E_0N}\mathcal{Z},$$ where $E_0$ is the ground state energy per spin, and $N$ is the number of spins. For the ferromagnetic case, this limit can be taken by absorbing the prefactor $\e^{\beta J}$ into the matrices $t$, and taking the limit on the level of the individual tensors. For the antiferromagnetic case, as a consequence of frustration, the same procedure yields matrices $t$ with elements that diverge as $\beta\to\infty$, so that the limit cannot be taken on the level of the local tensors and we cannot probe the $T=0$ partition function directly. Moreover, this inability to ‘regularize’ the partition function on the level of the local tensors makes the contraction of the tensor network problematic even at low temperatures. This is illustrated in Fig. \[fig:vumps\] where we show the performance of standard contraction methods for the frustrated Ising model ($J>0$) at different temperatures.

In this paper we address this issue by constructing a generic framework to map frustrated models defined on a trivial phase space to frustration-free models on a correlated phase space. This will enable us to understand the ground state ensemble as characterised by *local rules*, and yields a natural tensor-network construction allowing us to study the $T=0$ and low-temperature physics.

![Convergence of vumps algorithm [@ZaunerStauber2018; @Fishman2018; @Vanderstraeten2019] for the antiferromagnetic triangular lattice Ising model at different inverse temperatures $\beta$ with the standard tensor-network construction (dots) and the minimally-frustrated one (crosses), each for MPS with a bond dimension of $\chi=80$. We use a variational convergence measure, see Ref. . We observed similar behaviour using the corner transfer matrix renormalization group algorithm [@Orus2009], and a similar issue was observed for real-space renormalization techniques in Ref. .[]{data-label="fig:vumps"}](vumps.pdf){width="\columnwidth"}

***Minimal frustration.*** We consider a system of $d$-level spins $s_i$ on a lattice, and a translation invariant Hamiltonian $\mathcal{H}$ with local interaction terms $h_n$ of strictly bounded range. Starting from some cluster of spins $u$, the central idea relies on covering the lattice with the set $\mathcal{T}_u$ of overlapping translations of $u$, as illustrated in Fig. \[fig:cluster\], such that the Hamiltonian can be rewritten as a sum of strictly local terms, all acting within a single cluster of $\mathcal{T}_u$. Such a *tessellated Hamiltonian* is not uniquely defined since some interactions may be shared between overlapping clusters. We are therefore free to associate to each shared Hamiltonian term a weight $\alpha_n$, specifying how much of the term is accounted for in each cluster. $$\begin{aligned}
\mathcal{H} &= \sum_{c \in \mathcal{T}_u} \sum_{n \in c} \,\alpha^{c}_{n} \,h_n = \sum_{c \in \mathcal{T}_u} H^{\{\alpha\}}_c \label{eq:hamtess}\\
& \sum_{c \in \mathcal{T}_u| n \in c} \alpha^{c}_n = 1 \,,\qquad\forall \,n
\label{eq:alpha}\end{aligned}$$

The $\alpha^{c}_{n}$ are chosen such that the $H^{\{\alpha\}}_u$ are translation invariant. In this form, the Hamiltonian contains only terms that act locally within clusters, and the energy of any given state will be the sum of the local energies. Thus, the minimum of $H^{\{\alpha\}}_u$ –with respect to all spin configurations $C$ of $u$– immediately implies a lower bound on the global ground state energy. This lower bound can be maximized as a function of the $\alpha^{u}_{n}$ [@Huang2016]. Moreover, if there is a state whose energy is equal to this lower bound, it must minimize each one of the local Hamiltonians, and be a ground state. Effectively, the model is then “frustration free”: taking the spin configurations on $u$ as our local degree of freedom, each Hamiltonian term $H^{\{\alpha\}}_u$ can be minimized by selecting the right local configurations. The model’s inherent complexity has been moved to the constraints coming from the overlapping clusters, making this a highly correlated phase space. If the energy lower bound is saturated, we will say that the Hamiltonian tessellation in Eq.  has *minimal frustration*.

![A tessellation of spins on the kagome lattice: the cluster $u$ consists of twelve spins, and it shares five spins with each of the translated clusters $T_x(u)$ and $T_y(u)$.[]{data-label="fig:cluster"}](diagrams.pdf){width="0.5\columnwidth"}

One can draw a useful analogy between the correlated phase space and tiling problems. For a given local Hamiltonian $H^{\{\alpha\}}_u$, one can consider the local configurations of the cluster $u$ as *tiles*, with the overlapping spins determining whether they fit. One can thus try to prove minimal frustration either by determining the ground state energy and checking that it is equal to the lower bound, or by tiling the plane using only ’tiles’ that minimise $H^{\{\alpha\}}_u$ locally. Both these approaches lead to famously undecidable problems [@Wang1961; @Robinson1971; @Wang1975; @Wang2017]. If a Hamiltonian tessellation has minimal frustration, *all* the ground states can be constructed by tiling the minimal energy configurations of $H^{\{\alpha\}}_u$, and determining the macroscopic degeneracy of the ground state boils down to counting tilings.

***Tensor network construction.*** The latter can be accomplished efficiently with tensor networks. Indeed, if we have minimal frustration, we can represent the partition function $\mathcal{Z}$ in close analogy to the standard construction in Eq. . The shape of the tensor network is such that every vertex coincides with a cluster of $\mathcal{T}_u$, where we place a delta tensor, representing the configurations of $u$, on each vertex. Next, we introduce matrices on the bonds between vertices to enforce that spins shared between two clusters are in the same state in both clusters, essentially encoding the correlated nature of the phase space. Finally, since we only have on-site interactions in the tessellated Hamiltonian, we can include the Boltzmann weights in the non-zero elements of the vertex tensors. The resulting tensor network has a bond dimension equal to the number of configurations of the local cluster. This can be decreased significantly by performing a singular-value decomposition on the rank-deficient bond matrices. It is now possible to take the limit $\beta\to\infty$, as this will simply put a number of Boltzmann weights to zero, thereby effectively removing all configurations of the local cluster whose local energy is not minimal— a step that also shrinks the bond dimension. This new tensor network allows for a direct study of the system at $T=0$.

***Preliminary example.*** In order to illustrate the above construction we now reconsider the triangular lattice Ising antiferromagnet, for which the standard tensor-network approach was earlier shown to fail. We choose $\mathcal{T}_u$ to be the set of all clusters that form a square: $$\diagram{2}.$$ Here, each cluster has four interaction terms shared with neighbouring clusters. Their associated weights are constrained to be $\alpha_{1,2}'=1-\alpha_{1,2}$ by Eq. . If we choose $\alpha_1=\alpha_2=1/2$, we find that $H^{\{\alpha\}}_u\geq-J$. The ground state energy per cluster is known to be $-J$, so this choice of weights lead to a Hamiltonian tessellation with minimal frustration. The configurations of $u$ that minimize $H^{\{\alpha\}}_u$ (i.e., the tiles) are those with only one pair of aligned spins per triangle, exactly what one would expect. If we now construct the tensor network with these tiles for inverse temperatures $\beta=2$ and $\beta=\infty$, the standard contraction techniques nicely converge (Fig. \[fig:vumps\]).

***Spurious tiles.*** It can happen that *some* of the tiles minimizing $H^{\{\alpha\}}_u$ cannot fit into any tiling. We call such tiles *spurious*. In addition to obscuring the ground state rules and making it harder to understand the $T=0$ ensemble, spurious tiles can even revive the initial convergence problem. For instance, this happens in the triangular case when we consider $\alpha_1=\alpha_2=1$, so that only the down-triangles will contribute to the energy. The lowest energy of $H^{\{\alpha\}}_u$ is still $-J$, so there is still minimal frustration, but now the up-triangle is unrestricted and thus we have a larger number of tiles. Necessarily, these additional tiles are spurious, simply because there are $\alpha_1$ and $\alpha_2$ for which they are not ground state tiles. But in general, determining whether some of the tiles are spurious can be extremely challenging.

***Implementation.*** As stated before, any Hamiltonian tessellation gives rise to a lower bound on the ground state energy. If we maximize this bound, with respect to $\alpha^{u}_{n}$, we get an optimal lower bound $E_u$ for the cluster $u$: $$E_u\quad\leftarrow\quad\max_{\vec{\alpha}} \,E\,, \quad \text{with} \quad \left\{ \begin{array}{l}
H^{\{\alpha\}}_u(C)\geq E \; \forall C \\ \sum_{c \in \mathcal{T}_u| n \in c} \alpha^{c}_n = 1 \end{array} \right.
\label{linprogE}$$ This problem has the form of a linear program [@Grotschel1988], and may be solved using a standard linear programming toolbox. There are however many solutions of \[linprogE\], namely all $\alpha^{u}_{n}$ satisfying $$H^{\{\alpha\}}_u(C)\geq E_u, \quad \text{for all configurations $C$ of $u$},
\label{eq:ineq}$$ In the space of $\alpha^{u}_{n}$, the set of $\alpha^{u}_{n}$ for which all these inequalities are satisfied takes the form of a convex set, which we will refer to as $A^u$. The set of ground state tiles is fixed in the interior of $A^u$, but on a point of the boundary there will appear extra minimal energy configurations since inequalities that define this boundary become equalities. These extra configurations must thus correspond to spurious tiles. To have a minimal number of such spurious tiles, we therefore need to find an $\alpha^{u}_{n}$ in the interior of $A^u$. Note though, that there could yet be spurious tiles for such an $\alpha^{u}_{n}$, but at least they will be minimal in number, for a given $u$.

To find a point in the interior of $A^u$, we first find a set of points on the boundary, that together form a simplex of the same effective dimension as $A^u$. We can then take a point in the interior of this simplex and be guaranteed this point also lies in the interior of $A^u$. The way to build such a simplex is described by Alg. \[alg:simplex\].

However, direct use of Eq.  and Alg. \[alg:simplex\] can quickly become intractable since the number of constraints scales exponentially in the number of spins per cluster. But, technically, of all the inequalities in Eq. , we only need the ones that define the corners of the interior simplex. We could thus imagine systematically and progressively incorporating inequalities as we build this simplex. The concrete steps one needs to take are depicted Alg. \[alg:lean\]. This algorithm provides the energy lower bound, as \[linprogE\] did, and constructs an interior simplex while using only a manageable number of inequalities.

The bottleneck of this method is a step where the lowest energy configurations for given $\alpha^{u}_{n}$ are needed. This can be done by simple iteration, mixed integer linear programming [@Huang2016; @Larsen2018], or even Monte Carlo, but there is no escaping that finding all those configurations is an NP-complete problem. Therefore, although the memory cost remains low, the computation time grows exponentially as larger clusters are considered. Note that there exist models that can’t have minimal frustration for any finite cluster. One such model is discussed in [@Wang2017]. It is an open question whether such models could ever exhibit residual entropy.

$\vec{R}\gets$ random vector $\vec{\alpha}_1$$\gets$ $\max \,\vec{R} \cdot\vec{\alpha}\quad$ with $\alpha\in A$ $\vec{\beta}$ $\gets$ a point in simplex$[\{\vec{\alpha}_i\}]$ Translate $\vec{\alpha}$-space by $-\vec{\beta}$ $\{\vec{w}_i\}$ $\gets$ a basis of orthogonal vectors to all $\{\vec{\alpha}_i\}$ $\vec{\alpha}$$\gets$ $\max \,\vec{v} \cdot\vec{\alpha}\quad$ with $\alpha\in A$ Add $\vec{\alpha}$ to the set $\{\vec{\alpha}_i\}$ Return to the top of the while loop Stop the while loop $\{\vec{\alpha}_i\}$

$\{c_i\}$ $\gets$ choose some random configurations\
*Add random configurations to $\{c_i\}$ until there\
is a finite $E$ and a finite interior simplex* Add $c$ to $\{c_i\}$ $E_{\text{temp}}\gets$ solve Eq. \[linprogE\] for configurations $\{c_i\}$ $\{\vec{\alpha}_i\}\gets$ Update interior simplex for $\{c_i\}$ Return to the top of the while loop Stop the while loop $E_{\text{temp}}$, $\{\vec{\alpha}_i\}$

***Kagome lattice.*** As a challenging test case we consider a frustrated Ising model $$\mathcal{H} = J_1\sum_{\langle ij \rangle}\sigma_i\sigma_j + J_2\sum_{\langle\langle ij \rangle\rangle}\sigma_i\sigma_j + J_3 \sum_{\langle\langle\langle ij \rangle\rangle\rangle}\sigma_i\sigma_j\,,
\label{hamiltonian}$$ on the kagome lattice, where the sums run over (distance-based) first, second, and third nearest neighbours respectively. We take $J_1=-1$, $J_2=J_3=10$, which is unlikely to represent any realistic system, but exhibits interesting zero-temperature features. As a base cluster $u$, we use a full kagome star (12 spins, Fig. \[fig:cluster\]). The aforementioned linear program provides us with a ground state energy lower bound as well as candidate ground state tiles (Fig. \[fig:tile\]). It is known that the ground state energy per site of this model is $E = \frac{2}{3} J_1 - \frac{2}{3} J_2 - J_3$, matching our result. We thus have a tessellated Hamiltonian with minimal frustration.

![The residual entropy per site of the minimally frustrated tensor network for the model in Eq. , obtained with the vumps algorithm with different bond dimensions $\chi$. In the inset we show that the value converges to the value for the Ising antiferromagnet on the triangular lattice [@Wannier1973].[]{data-label="fig:entropy"}](entropy.pdf){width="0.99\columnwidth"}

The ground state tiles that we get can be classified in two types according to the spin configuration on the hexagon: 4 up/down-spins and 2 down/up spins for type-I, same number of up-spins and down-spins for type-II, see Fig. \[fig:tile\]. We proceed to understanding the ground state ensemble of this model by first characterizing the type-I ensemble, and then describing how type-II tiles modify this picture.

The ensemble generated by type-I tiles is best understood in terms of a mapping to a loop model, where we imagine a line between anti-aligned spins, as indicated in Fig. \[fig:tile\]. Clearly the lines cannot merge or end in a tile, so the lines must form loops. In any given type-I tile, either all the up triangles have a line run through them, or all the down triangles do. The way tiles have to fit together makes it clear that the lines can only go through one kind of triangle in any one configuration, so there must be reflection-symmetry breaking. The lines live on a honeycomb lattice which is generated by connecting the centre of each hexagon to its neighbouring up triangles (respectively down triangles), and each site of this honeycomb lattice has to be crossed by a loop. We thus end up with a loop model that is the complement to the dimer model on the honeycomb lattice, where the dimers would indicate the absence of a line. The latter is itself dual to the triangular-lattice Ising antiferromagnet, for which an exact value for the residual entropy is known [@Wannier1973].

The value for the entropy that we find numerically (see inset of Fig. \[fig:entropy\]) suggests that the type-II tiles are not relevant for the entropy in the thermodynamic limit. In finite-size systems, however, type-II tiles show up as domain walls between different sectors of reflection symmetry (Fig. \[fig:domainwall\]). They are thus not spurious, but they lead to a sub-extensive contribution to the residual entropy. This point will be discussed further in a future work devoted to the general $J_1$-$J_2$-$J_3$ model on the kagome lattice [@preparation].

![An example of a spin configuration in the ground state, and the corresponding red line configuration. Type-II tiles form domain walls between the symmetry broken sectors formed by type-I tiles. This configuration was generated by Monte Carlo on a finite system with periodic boundary conditions using a combination of the cluster update from Ref. , single spin flip and parallel tempering.[]{data-label="fig:domainwall"}](domain_wall.pdf){width="0.7\columnwidth"}

***Outlook.*** In this paper, we have provided a generic framework for studying frustrated spin systems, identifying effective degrees of freedom in these models at low temperature. The concept of minimal frustration, both Alg. \[alg:simplex\] and Alg. \[alg:lean\] as well as the tensor network construction work in all dimensions. It might be possible to devise a Monte Carlo method that leverages the tools developed in this work, perhaps in combination with Ref. .

The issue of contracting a tensor network with very large and very small numbers is strongly reminiscent of the infamous sign problem in Monte Carlo. It seems that this issue can pop up in any tensor network, in particular in PEPS. Here we solve an avatar of this “sign problem for tensor networks”, potentially showing the way to a more general solution.

Finally, we can consider the effect of quantum dynamics on the correlated phase spaces. Indeed we can write down PEPS wavefunctions by promoting the tiles to quantum degrees of freedom to effectively describe quantum corrections that would be present in any real life material.

***Acknowledgements.*** This work was initiated at the Centro de Ciencias de Benasque Pedro Pascual. FV thanks Michael Lawler for discussions about tensor networks and frustrated systems. This research is supported by ERC grant QUTE (647905) and FWO (G0E1820N) (LV, BV, FV), and the Swiss National Science Foundation (JC, FM).
