# Critical line of the triangular Ising antiferromagnet in a field from a
$C_{3}$-symmetric corner transfer matrix algorithm

Samuel Nyckees<EMAIL_ADDRESS>Institute of Physics, Ecole
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Afonso Rufino  Institute of Physics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Lausanne, Switzerland  Frédéric Mila  Institute of Physics,
Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Jeanne Colbois<EMAIL_ADDRESS>Laboratoire de Physique Théorique,
Université de Toulouse, CNRS, UPS, 31400 Toulouse, France

(August 28, 2024)

###### Abstract

The corner transfer matrix renormalization group (CTMRG) algorithm has been
extensively used to investigate both classical and quantum two-dimensional
(2D) lattice models. The convergence of the algorithm can strongly vary from
model to model depending on the underlying geometry and symmetries, and the
presence of algebraic correlations. An important factor in the convergence of
the algorithm is the lattice symmetry, which can be broken due to the
necessity of mapping the problem onto the square lattice. We propose a variant
of the CTMRG algorithm, designed for models with $C_{3}$-symmetry, which we
apply to the conceptually simple yet numerically challenging problem of the
triangular lattice Ising antiferromagnet in a field, at zero and low
temperatures. We study how the finite-temperature three-state Potts critical
line in this model approaches the ground-state Kosterlitz-Thouless transition
driven by a reduced field ($h/T$). In this particular instance, we show that
the $C_{3}$-symmetric CTMRG leads to much more precise results than both
existing results from exact diagonalization of transfer matrices and Monte
Carlo.

##  I Introduction

A very successful approach to study classical spin models on two dimensional
lattices is to write the partition function of the system as the contraction
of a tensor network [1, 2, 3, 4, 5]. Inspired by the original transfer matrix
formulation [6, 7], this tensor network associates to each interaction a
tensor carrying its Boltzmann weight and is typically defined on the same
lattice as the original model. With the exception of a few well-known exactly
solvable cases, namely zero-field Ising models on planar lattices [7, 8, 9,
10, 11, 12], an exact evaluation of this partition function is generally
exponentially hard. In those cases, an approximate contraction scheme, such as
tensor network renormalization group (TRG, TNR) [4, 13], boundary matrix
product state (MPS) with infinite time-evolved block decimation (iTEBD) [14,
15, 16], variational uniform MPS (VUMPS) [17, 18, 19], or the corner transfer
matrix renormalization group (CTMRG) [1, 2, 5, 18] can provide very accurate
results.

With the notable exception of TRG/TNR [4, 20, 21, 22], these algorithms are,
in most cases, formulated on the square lattice. In particular, this is the
case for CTMRG, which has become a cornerstone of tensor network approaches in
classical and quantum lattice models. This algorithm was first introduced by
Nishino and Okunishi [1, 2, 3] as an efficient contraction scheme to evaluate
partition functions of infinite two-dimensional square lattice models. It
comes as a combination of Baxter’s corner transfer matrix [23, 24, 25, 26] and
White’s density matrix renormalization group (DMRG) algorithm [27, 28] .
Although initially introduced in the context of 2D statistical physics as a
complementary approach to other transfer matrix methods or to Monte Carlo, its
use rapidly extended to 2D quantum physics. Indeed, in the last two decades,
CTMRG has been used as a contraction algorithm for infinite projected
entangled pair states (iPEPS) wavefunctions [29, 15, 5, 30, 22]. Pairing CTMRG
with different update schemes for iPEPS such as simple [31] and full update
[32] and more recently automatic differentiation [33, 34] has notably
supported new results in fermionic systems [35, 36, 37, 38] and frustrated
systems [39, 40, 41, 42].

Besides the now well-established power of tensor networks for quantum systems,
a promise of tensor networks for classical spin systems has been to offer
either a powerful alternative to Monte Carlo (see e.g. Refs. 43, 44, 45, 46,
47) or a support to improve sampling [48, 49]. However, in recent years, their
application to frustrated two- and three-dimensional classical models has
attracted some attention [50, 51, 52, 53, 54, 55, 56, 57, 58, 59], as it has
been established that particular care has to be taken in the tensor network
formulation to avoid numerical instabilities or convergence to wrong results.
Importantly, these instabilities seem to get significantly reduced away from a
macroscopically degenerate ground state [50], but they can create issues close
to a critical ground state such as that of the triangular lattice Ising
antiferromagnet [53]. In the context of exact, differentiable contraction,
these instabilities can be dealt with by working with the logarithm of the
Boltzmann weight [54]. In the case of approximate contraction, where it is
unclear how to use such a construction, a successful alternative approach has
been to rely on ensuring that the ground-state local rule is satisfied at the
level of the tensor, such that the low-temperature limit of the tensor network
remains well-defined [53, 55, 58]. This approach has been mostly validated
using VUMPS; here, we find that it is also successful with various CTMRG
algorithms.

Since CTMRG is defined on the square lattice, a usual first step is to map the
tensor network from the original lattice onto the square lattice [60, 61, 62,
63, 64, 65, 42]. While it has been generally successful (of particular
interest here, there is the case of the triangular lattice Ising ferromagnet,
[60]), this step may occasionally lead to problems such as poor convergence or
non-physical breaking of symmetry [66]; in contrast, it is well known that
improved performance is achieved when the algorithm makes use of the
underlying symmetry of the lattice [67, 66]. In this spirit, proposing
variants of existing algorithms formulated to fit the lattice symmetry of the
problem could allow investigation of systems with better accuracy.

In this paper, we design a CTMRG algorithm to contract infinite tensor
networks defined on the honeycomb lattice, naturally giving rise to a
$C_{3}$-symmetric contraction scheme. A simple construction enables us to
apply it to the low temperature phase diagram of the triangular
antiferromagnet Ising model in a field. In Section II we describe this model
and recall previous works on the topic. In Section III, we recall the CTMRG
algorithm on the square lattice and then introduce the CTMRG algorithm on the
honeycomb lattice. In Section IV, we first benchmark the algorithm on the
classical Ising antiferromagnet on the triangular lattice, naturally
formulated through a dual construction as a tensor network on the honeycomb
lattice. We then revisit the effect of the magnetic field both in the
constrained model and the finite-temperature case. We show that the location
of the critical line at finite temperature can be evaluated with much higher
precision than previously achievable. Finally, in Section V we discuss our
results and provide an outlook.

##  II The model

The antiferromagnetic Ising model on the triangular lattice in a magnetic
field is defined by the Hamiltonian:

$\mathcal{H}=J\sum_{\langle
i,j\rangle}\sigma_{i}\sigma_{j}-h\sum_{i}\sigma_{i}$ (1)

with $\sigma\in\\{+1,-1\\}$ an Ising spin variable, $J>0$ the interaction
parameter and $h$ the field parameter.

We show a sketch of its phase diagram in Fig. 1. In the absence of a magnetic
field, the model has a macroscopically degenerate ground-state [10, 11], i.e.
a finite residual entropy, and is characterized by each triangle satisfying a
two-up one-down, two-down one-up (UUD/DDU) rule for the spins. This critical
point has algebraically decaying spin-spin correlation [68] characterised by a
critical exponent $\eta=1/2$ and central charge $c=1$, and is commonly
referred to as the Villain-Stephenson (VS) point. It can be described using a
Coulomb-gas construction [69, 70, 71] thanks to an exact mapping onto a
triangular solid-on-solid (SOS) model [72, 70, 73]. Upon introducing a
positive field, the configurations with each triangle having UUD spins are
favored, giving rise to long-range correlations. This $m=1/3$ magnetization
plateau, with $\sqrt{3}\times\sqrt{3}$ symmetry breaking stabilizes until
$h=6J$. This point has a finite residual entropy $S=0.333242...$ [74, 75]. For
$h>6J$, the ground state is fully polarized.

Figure 1: (a) Sketch of the phase diagram in the $(h/J,T/J)$ units. The
ordered phase melts through a three-state Potts transition. The slope of the
critical line near the $m=1$ magnetization plateau is determined by a mapping
to the hard hexagon model. (b) Phase diagram of the model in reduced
coordinates $(H,K^{-1})$ with $H=h/T$ and $K=J/T$. At $K^{-1}=0$, the system
maps onto an SOS model with a critical phase and an ordered phase separated by
a KT transition.

At zero field, for any nonzero temperature, the system has a finite
correlation length which diverges exponentially fast as $T\rightarrow 0$ [76,
77]. In contrast, for $0<h<6J$, the $\sqrt{3}\times\sqrt{3}$ ordered phase
melts at finite temperature through a three-state Potts transition as expected
from the ground-state symmetry [78] and verified by phenomenological scaling
and transfer matrix calculations [79, 80, 81, 82, 83, 84]. This transition is
characterized by critical exponents $\beta=1/9,\nu=5/6,\eta=4/15,c=4/5$, with
$\beta$ describing the order parameter power law in the ordered phase and
$\nu$ characterizing the divergence of the correlation length.

The shape of the Potts critical line in the limit $T\rightarrow 0$ has
attracted quite some interest. In the vicinity of $h=6J$, the shape of the
line is predicted by the hard-hexagon model whose second-order phase
transition is also in the three-state Potts universality class. Indeed, the
partition function of the triangular Ising antiferromagnet in the limit
$T\rightarrow 0,h\rightarrow 6J$ with $(6J-h)/T=:\ln(\zeta)/2$ remaining
finite maps to that of the hard hexagon model with fugacity $\zeta$ [85, 84];
thus the slope of the three-state Potts critical line in that limit is
determined by the critical fugacity of the hard-hexagon model
$\zeta_{c}=\frac{11+5\sqrt{5}}{2}$ [86] as

$T_{c}\rightarrow\frac{1}{\ln(\zeta_{c})}(12J-2h)$ (2)

The limit $T\rightarrow 0,h\rightarrow 0$ has proven much more challenging. A
first conjecture [79] was that the Potts line would approach the VS point with
an infinite slope, but renormalization group investigation instead showed that
in this limit, the slope of the transition line should remain finite [70]. The
model is sometimes introduced considering the reduced Hamiltonian with
parameters $H=h/T$ and $K=J/T$:

$\mathcal{H}_{R}=\frac{\mathcal{H}}{T}=K\sum_{\langle
i,j\rangle}\sigma_{i}\sigma_{j}-H\sum_{i}\sigma_{i}$ (3)

where an additional critical phase appears at $K^{-1}=0$. Indeed, introducing
a finite reduced field at the VS point acts as a perturbation which become
relevant only when $\eta=4/9$ where a Kosterlitz-Thouless (KT)[87, 88, 70]
transition into an ordered phase must occur. The location of the KT transition
is not predicted by the RG analysis, but transfer matrix studies have located
it around $H_{KT}=0.266\pm 0.01$ [82, 83, 84, 89, 90]. The three-state Potts
transition should then meet the KT point in the zero temperature limit. A
study based on transfer matrix and renormalization group analysis [84] has
given evidence in favor of that scenario and further suggested that the
transition line in the $(K^{-1},H)$ phase diagram approaches $H_{KT}$ with a
square root singularity. However the smallest accessible critical field of the
three-state Potts transition remains significantly larger than $H_{KT}$ and a
definitive answer on the shape of the transition near $H_{KT}$ remains to be
given.

For the rest of the paper, we will sometimes refer to $K^{-1}$ as the
temperature but will always distinguish between the reduced field $H$ and the
field $h$.

##  III Corner transfer matrix renormalization algorithm

In this section we give a brief overview of the square CTMRG and its
implementation for the triangular Ising antiferromagnet. We then introduce a
new version of the CTMRG algorithm that contracts honeycomb tensor networks
and show how to express the partition function of the antiferromagnetic Ising
model on the triangular lattice as the contraction of an infinite honeycomb
tensor network.

###  III.1 Algorithm for the square lattice

The square CTMRG algorithm approximates the contraction of infinite square
tensor networks made of rank-four local tensor $a$ of dimension $d\times
d\times d\times d$ with a contraction of eight different tensors referred to
as the environment and the local tensor $a$ as shown in Fig. 2. The
environment is made of four corner tensors $C_{i}$ of dimension
$\chi\times\chi$ and four edge tensors $T_{i}$ of dimension $\chi\times
d\times\chi$. The approximation is controlled by the bond dimension $\chi$ and
in the infinite bond dimension limit one recovers the exact result.

Figure 2: The partition function written as the infinite contraction of a rank
four tensor in a square lattice is approximated by a contraction of nine
different tensors. Bold lines represent legs of bond dimension $\chi$ while
thin lines represent bond dimension $d$.

In order to find a representation of the infinite lattice tensor network with
a finite bond dimension, tensors are iteratively added into the network and
the thermodynamic limit is obtained when some observables have converged.
There are different ways of constructing the environment and in the present
case we focus on the directional CTMRG algorithm [5], where columns are added
to the network one at a time. It consists of two steps, which when repeated
along the four lattice directions define a full CTMRG iteration:

  * •

Absorption: a new column is introduced into the network in one of the
directions and contracted with the appropriate corner and edge tensors,
increasing their dimension by a factor of $d$. An illustration of the
absorption for a left move is included in Fig. 3.

  * •

Renormalization: if one of the dimensions of a tensor exceeds the cutoff
dimension, the dimension of this index is truncated after application of an
isometry $\mathcal{U}$. This step is illustrated in Fig. 4 for a left move.

Figure 3: Illustration of the absorption step of square CTMRG in tensor
network notation. Figure 4: Illustration of the renormalization step of square
CTMRG in tensor network notation.

There are several possible choices for the definition of the isometry
$\mathcal{U}$, which can affect the convergence of the algorithm. In the
present case we used the isometry first introduced by Corboz in [91]. 111The
isometries first introduced by Nishino and Okunishi in Ref. [1] or Orus in
Ref. [5] performed poorly at small reduced field and close to the critical
regime.

One of the advantages of the CTMRG formalism is that it gives access to the
entanglement entropy, whose divergence at second order phase transitions can
be used to study the associated universality classes [92, 44, 93] and is
defined as :

$\displaystyle S$ $\displaystyle=-\text{Tr}(\rho\log(\rho))$ (4)

by setting the density matrix $\rho$ as:

$\displaystyle\rho$
$\displaystyle=C_{1}C_{2}C_{3}C_{4}/\text{Tr}(C_{1}C_{2}C_{3}C_{4}).$ (5)

#### Mapping on the triangular lattice

For a given statistical mechanical system the choice of a tensor network
representation of its partition function is not unique. In the case of
frustrated systems such as the triangular Ising antiferromagnet the choice of
the tensor greatly affects the convergence. It was observed [53, 56] that it
is best to use a tensor network construction obtained by factorizing the
Boltzmann weight of a global spin configuration into a product over a set of
tiles whose reunion covers the whole lattice. In this spirit, we define a
tensor for the triangular Ising antiferromagnet that accounts for the
Boltzmann weight of a hexagon of seven spins given by :

$\displaystyle
a_{(\sigma_{1},\sigma_{2}),(\sigma_{7},\sigma_{3}),(\sigma_{5},\sigma_{4}),(\sigma_{6},\sigma_{8})}=\sum_{\sigma_{0}=\pm
1}e^{-\beta
H(\sigma_{0},\dots,\sigma_{6})}\delta_{\sigma_{2},\sigma_{7}}\delta_{\sigma_{5},\sigma_{8}}$
(6)

with

$H(\sigma_{0},\dots,\sigma_{6})=J\left(\sigma_{0}\sum_{i=1}^{6}\sigma_{i}+\sum_{i=1}^{5}\sigma_{i}\sigma_{i+1}+\sigma_{6}\sigma_{1}\right)\\\
-h\sigma_{0}-\frac{h}{3}\sum_{i=1}^{6}\sigma_{i}$ (7)

We give the schematic picture of the tensor in Fig. 5.

Figure 5: Local tensor made of seven spins given by Eqs. (6) and (7).

###  III.2 Algorithm for the honeycomb lattice

We now introduce a variant of CTMRG that contracts infinite tensor networks
defined on the (bipartite) honeycomb lattice with local tensors $a$ and $b$ of
dimension $d\times d\times d$. It is worth noting that there are other CTMRG
algorithms on lattices other than square [94, 95]. The most relevant in our
case has originally been introduced by Gendiar et al. [96] while studying
Ising models on triangular-tiled hyperbolic lattices formulated as
interactions-round-a-face (IRF) tensor networks. This was recently adapted to
the context of 2D-quantum physics model using automatic differentiation [66]
with iPEPS wave-functions defined on a honeycomb lattice with a single local
tensor $a$. Although they apply to similar lattices, we see a key difference
between this algorithm and the one presented here. In Ref. 66, the corner
matrix represents a sixth of the partition function while in ours it
corresponds to a third. This in turn leads to inserting the isometries in
different places and to a more natural generalization of our algorithm when
the bipartite lattice is made of two different local tensors. We postpone this
discussion to the end of the present section.

Figure 6: Partition function written as a two dimensional honeycomb tensor
network, and then contracted to an environment made of nine tensors, and a
unit cell made of six on-site tensors. Bold lines represents bond dimension of
$\chi$ while thinner lines indicate a smaller bond dimension of $d$.

The algorithm presented here approximates the infinite tensor network with
nine tensors which represent the environment $E=\\{T_{\alpha}^{i},C_{i}\mid
i\in\\{1,2,3\\},\alpha\in\\{a,b\\}\\}$ as shown in Fig. 6 where the corner
tensors $C_{i}$ are of dimension $\chi\times\chi$, and the edge tensors
$T_{\alpha}^{i}$ are of dimension $\chi\times d\times\chi$. Just as in the
square lattice CTMRG case [2], the approximation is controlled by the bond
dimension $\chi$ and in the infinite $\chi$ limit, one recovers the exact
result. The algorithm iterates two alternating steps as well, the update and
the renormalization, until the convergence of some observable of
interest222See e.g. Refs. 22, 97 for a discussion of finite-size vs finite-
bond dimension in the square lattice CTMRG.. For the CTMRG on the honeycomb
lattice, we propose the following update (the full iteration is shown in Figs.
7 and 8):

$\displaystyle C_{i}^{\prime}$ $\displaystyle=abC_{i}T_{a}^{i}T_{b}^{f(i-1)}$
(8) $\displaystyle T_{b}^{i^{\prime}}$ $\displaystyle=T_{a}^{i}b$ (9)
$\displaystyle T_{a}^{i^{\prime}}$ $\displaystyle=aT_{b}^{i}$ (10)

where we have introduced $f(i)=\text{mod}(i,3)$ to make the notation lighter.
We then have to project the tensors in a relevant subspace in order to keep
their dimensions under control. This is called the renormalization step and
goes as:

$\displaystyle C_{i}$
$\displaystyle=\mathcal{U}_{i}C_{i}^{\prime}\mathcal{U}_{f(i-1)}^{\dagger}$
(11) $\displaystyle T_{b}^{i}$
$\displaystyle=T_{b}^{i^{\prime}}\mathcal{U}_{i}$ (12) $\displaystyle
T_{a}^{i}$ $\displaystyle=\mathcal{U}_{i}^{\dagger}T_{a}^{i^{\prime}}$ (13)

where the isometry $\mathcal{U}_{i}$ is computed by truncating the singular
value decomposition of $C_{f(i+1)}C_{f(i+2)}C_{i}$. This isometry would be the
equivalent in the square CTMRG of the isometry originally used by Nishino and
Okunishi [2]. We illustrate picturally in Fig. 9 how the iterations correspond
to contracting the partition function on the whole honeycomb lattice.

Figure 7: One full iteration for the corner tensor $C_{1}$, Eqs. (8) and (11).
$C_{2}$ and $C_{3}$ are updated similarly Figure 8: One full iteration for the
edge tensors $T_{a}^{1}$ and $T_{b}^{1}$, Eqs. (9, 10, 12, 13). The other edge
tensors are updated similarly. Figure 9: Schematic way of representing how the
environment is updated. Primed tensors represent the tensors after one
iteration.

Similarly to the square CTMRG we can compute the entanglement entropy by
setting the density matrix $\rho$ as:

$\displaystyle\rho$
$\displaystyle=C_{1}C_{2}C_{3}/\text{Tr}(C_{1}C_{2}C_{3}).$ (14)

A particularly useful case is when the local tensors $a$ and $b$ are invariant
by rotation. Then, one can reduce the environment to only three different
tensors, $E=\\{C,T_{a},T_{b}\\}$ and the density matrix becomes:

$\displaystyle\rho$ $\displaystyle=C^{3}/\text{Tr}(C^{3})$ (15)

We will refer to the algorithm making use of that symmetry as the $C_{3}$
symmetric (honeycomb) CTMRG.

In the case of local tensors which are symmetric under rotation and reflection
(without any particular relation between $a$ and $b$), the $C_{3}$ symmetric
CTMRG algorithm gives direct access to the transfer matrix whose spectrum can
be used to compute the correlation lengths and wave-vector in their respective
directions. Indeed, by denoting the normalized eigenvalues of the transfer
matrix as:

$\displaystyle\lambda_{i}=e^{-\epsilon_{i}+i\phi_{i}}$ (16)

one can show that the correlation length and wave-vector are given by:

$\displaystyle\xi_{\chi}=\epsilon_{2}^{-1}$ $\displaystyle\qquad
q_{\chi}=\phi_{2}.$ (17)

It was suggested by Rams et al. [98] to extrapolate the finite bond dimension
correlation length with respect to higher gaps $\delta$ in the transfer
matrices as:

$\displaystyle\xi_{\chi}^{-1}=\xi^{-1}_{\infty}+a\delta$ (18)

Indeed, when the bond dimension increases, the spectrum above the first gap
needs to converge to a continuum and $\delta$ goes to zero. We will be using
in the next section $\delta=\epsilon_{4}-\epsilon_{2}$. In Fig. 10, we show
the transfer matrix in the $x$ direction with the convention of direction
shown in Fig. 6.

Furthermore, if the infinite corner tensor is hermitian, one can keep $C$
diagonal and use an eigenvalue decomposition for better numerical stability.
In particular, this is the case if the tensor network has a $C_{3v}$ symmetry.

Figure 10: Transfer matrix in the $x$ direction in the case of fully symmetric
tensors $a$ and $b$.

We now briefly discuss the practical differences between the algorithm
introduced in Ref. 66 and the $C_{3}$-symmetric CTMRG introduced here. As
already mentioned, the main difference between the two algorithms is the
definition of the corner transfer matrix which in turn leads to a different
truncation scheme. In the present case we consider the CTM to represent a
third of the partition function while in Refs. 66 it represents a sixth. Using
Gendiar’s algorithm on a bipartite lattice, one would then need to define two
different corners $C_{a}$ and $C_{b}$ which in turn at the renormalization
step would naively require two different isometries $U_{1}$ and $U_{2}$ to be
projected onto their relevant subspace such that
$C^{\prime}_{a}=U_{1}C_{a}U_{2}$ and $C^{\prime}_{b}=U_{2}C_{b}U_{1}$. In
contrast using only one corner matrix $C$ which effectively represents
$C=C_{a}T_{a}T_{b}C_{b}$ only requires one isometry and $C^{\prime}=UCU$.

#### Mapping from the triangular to the honeycomb lattice

To apply the $C_{3}$-symmetric CTMRG algorithm to the triangular lattice Ising
antiferromagnet we map the triangular lattice onto the honeycomb lattice by
defining the local tensors $a$ and $b$ on the dual as:

$\displaystyle
a_{(\sigma_{1},\sigma_{4})(\sigma_{2},\sigma_{5})(\sigma_{3},\sigma_{6})}$
$\displaystyle=\delta_{\sigma_{1}\sigma_{6}}\delta_{\sigma_{4}\sigma_{2}}\delta_{\sigma_{3}\sigma_{5}}e^{\frac{H}{6}(\sigma_{1}+\sigma_{2}+\sigma_{3})}$
$\displaystyle\times
e^{-\frac{K}{2}}e^{-\frac{K}{2}(\sigma_{1}\sigma_{2}+\sigma_{2}\sigma_{3}+\sigma_{3}\sigma_{1})}$
(19) $\displaystyle
b_{(\sigma_{1},\sigma_{4})(\sigma_{2},\sigma_{5})(\sigma_{3},\sigma_{6})}$
$\displaystyle=\delta_{\sigma_{1}\sigma_{5}}\delta_{\sigma_{2}\sigma_{6}}\delta_{\sigma_{4}\sigma_{3}}e^{\frac{H}{6}(\sigma_{1}+\sigma_{2}+\sigma_{3})}$
$\displaystyle\times
e^{-\frac{K}{2}}e^{-\frac{K}{2}(\sigma_{1}\sigma_{2}+\sigma_{2}\sigma_{3}+\sigma_{3}\sigma_{1})}$

where $\delta$ denotes the Kronecker delta with $K$ and $H$ the reduced
coupling constants. We give a diagrammatic expression of the tensors $a$ and
$b$ in Fig. 11. A more detailed explanation on how to express a classical
partition function as the contraction of a two-dimensional network and how to
map the problem on a specific lattice geometry can be found in [53]. Defined
as such, the local tensors $a$ and $b$ are rotation invariant, and we can use
the $C_{3}$ CTMRG algorithm. We note that the tensors $a$ and $b$ are related
by $a_{i,j,k}=b_{k,j,i}$. Although the tensors $a$ and $b$ are not fully
symmetric and the previous construction of the transfer matrix (Fig. 10)
cannot be applied, we can still construct a transfer matrix but in a different
way as shown in Fig. 12 with $\delta$ a $4\times 4$ matrix representing the
permutation of two legs given by:

$\displaystyle\delta=$ $\displaystyle\begin{pmatrix}1&0&0&0\\\ 0&0&1&0\\\
0&1&0&0\\\ 0&0&0&1\end{pmatrix}$ (20)

This construction applies to the case of the triangular Ising model and does
not generalize to arbitrary rotational invariant tensors $a$ and $b$.

Figure 11: Definition of the local tensors $a$ and $b$. By grouping the
indices pairwise we end up with tensors $a$ and $b$ of dimension $4\times
4\times 4$.  Figure 12: Transfer matrix in the $x$ direction for the Ising
model for a $6\times\infty$ system on the triangular lattice. The same
reasoning applies to an infinite system in both directions and give the same
results. The red lines indicate lines that must be contracted. The dashed
lines represent the perspective in 3D and indicate that the legs pass below
the full lines.

In the ordered phase, the system undergoes a translational symmetry breaking
and the triangular lattice is divided into three sub-lattices $A,B$ and $C$
with different magnetizations (Fig. 1). We note that the honeycomb CTMRG
naturally allows for this freedom while imposing the translational symmetry
associated with the long-range correlations in the ground state. This property
can be deduced by looking at Fig. 13, where we have highlighted the different
sub-lattices in the local tensors $a$ and $b$. Indeed, any observable measured
on a site of the sub-lattice $A,B$ or $C$ (represented by colored dots) has to
be the same if measured on a different site of the same sub-lattice due to the
network being invariant by $C_{3}$. This is only valid when the corner and
edge tensors are identical in all three possible directions. As we will see,
enforcing this symmetry will lead to a better convergence and will allow us to
access lower temperatures than previously available.

Figure 13: Illustration of the implication of the $C_{3}$ symmetry on the
expectation value of local observables. The tensors $a$ and $b$ are given by
Eq. (III.2), and each tensor describes three sites. The $C_{3}$ symmetry in
the environment imposes that different sites on the same sub-lattice have the
same local magnetization. The figure highlights the sub-lattices on the
tensors $a$ and $b$. One site belongs to the sub-lattice A ( ). Three sites
belong to the sub-lattice B ( ) and three more sites belong to the sub-lattice
C ( ).

##  IV Results

We now present the results obtained by the $C_{3}$-symmetric CTMRG. We
benchmark the algorithm on the antiferromagnetic triangular Ising model at
zero and finite temperature. We then turn to the constrained model, where we
investigate the location of the KT transition. Finally, we discuss and compare
results obtained from the square and honeycomb CTMRG. We define the order
parameter as:

$\displaystyle\psi=\frac{1}{3}|\psi_{A}+e^{2\pi i/3}\psi_{B}+e^{-2\pi
i/3}\psi_{C}|$ (21)

with $\psi_{i}$ the magnetization on the sub-lattice $i$. In the ordered
phase, due to the $\sqrt{3}\times\sqrt{3}$ unit cell, the honeycomb CTMRG
algorithm converges to three different environments $E_{A},E_{B}$ and $E_{C}$.
Observables such as the magnetization in the middle of the plaquette or the
entanglement entropy thus converge modulo three. Similarly, the transfer
matrices obtained from $E_{i}$’s do not have the same spectrum. One cannot
then consider this as the proper transfer matrix and would need to take into
account all environments in order to construct the true transfer matrix. To do
so, we would need to generalize the algorithm to a CTMRG with a general,
multi-site unit cell [5]. We thus only have access to the correlation length
in the disordered phase. Yet, as we will see, a two-site unit cell CTMRG is
powerful enough for our purposes.

###  IV.1 Benchmark at zero field

As mentioned in the introduction, at $h=0$ we recover the solvable triangular
lattice Ising antiferromagnet whose energy is exactly known [10, 99, 11]. Its
ground state has critical correlations, and upon decreasing the temperature,
the correlation length is known to diverge exponentially fast [76, 77] as :

$\displaystyle\frac{1}{\xi}=-\log(\tanh(\beta J))$ (22)

and this critical point belongs to the Villain-Stephenson universality class
[76, 100, 101]. We compare our results with the exact values in Fig. 14 and 15
where the energy and correlation length are in good agreement with the
associated theoretical predictions.

Figure 14: Simulations were done at $h=0$ for two different bond dimensions.
The error is of the order $10^{-13}$ for temperatures larger than 0.6 where
the correlation length is small. For the smallest computed temperature
($T=0.1$), the error is of the order $10^{-10}$. Figure 15: Simulation done at
$h=0$. Top panel : the correlation length with respect to the temperature. The
straight line represents a linear fit of the inverse correlation length.
Bottom panel : the extrapolation used for the correlation length with respect
to the gap $\delta=\epsilon_{4}-\epsilon_{2}$ in the transfer matrix. We have
used bond dimensions $\chi\in[100,200]$.

At zero field, the system is critical with infinite correlation length.
However, when introducing a finite bond dimension $\chi$, the correlation
length will become finite as well and can only be increased by considering a
larger bond dimension. It has been argued that it follows a power law
$\xi\sim\chi^{\kappa}$ with universal exponent $\kappa$ given by [102, 103]:

$\displaystyle\kappa=\frac{6}{c(\sqrt{12/c}+1)}.$ (23)

Using the correlation length scaling, one can derive the following finite bond
dimension dependency of the order parameter and entanglement entropy [102,
103] :

$\displaystyle\psi=\chi^{-\frac{\kappa\eta}{2}}$ (24) $\displaystyle
S=\frac{c\kappa}{6}\log(\chi).$ (25)

Assuming the central charge, we can thus define effective exponents $\kappa$
and $\eta$ by fitting the order parameter and entanglement entropy with Eqs.
(24) and (25) over a certain range of bond dimensions $\mathcal{I}_{\chi}$. By
defining $\alpha$ and $\beta$ such that

$\displaystyle S=\beta\log(\chi)$ (26)
$\displaystyle\log(\psi)=\alpha\log(\chi).$ (27)

and assuming the central charge we then get $\kappa=6\beta/c$ and
$\eta=-2\alpha/\kappa$ and the computation of the critical exponent is reduced
to a simple linear fit. As the environment converges modulo three, the
entanglement entropy takes three different values $S_{A},S_{B},S_{C}$
depending on the modulo of the number of iterations. Furthermore, by symmetry
of the model we have $S_{B}=S_{C}$ and we can then define two different
exponents $\eta_{A}\equiv\eta(S_{A})$ and $\eta_{B}\equiv\eta(S_{B})$
depending on which entropy $S_{A}$ or $S_{B}$ we choose to define $\kappa$
from. We benchmark the methodology at the Villain-Stephenson point at
$H=K^{-1}=0$ using bond dimension $\chi\in[200,400]$ and compare the results
with the square CTMRG. The results are shown in Fig. 16. Assuming $c=1$, we
found for the honeycomb CTMRG $\eta_{A}=0.4990\pm 0.0027$ and
$\eta_{B}=0.4986\pm 0.0027$ in good agreement with the exact value $\eta=1/2$
where the error bars have been computed from the goodness of the fit.
Conversely, one could have assumed $\eta=1/2$ and we recover $c_{A}=1.002\pm
0.005$ and $c_{B}=1.003\pm 0.005$. However, we note that even by assuming the
central charge $c=1$ we find $\kappa_{A}=1.230\pm 0.003$ and
$\kappa_{B}=1.231\pm 0.003$ which do not agree with the theoretical prediction
given by Eq. (23). We attribute this behavior to various corrections in bond
dimension [22, 104, 44]. We then choose to assume $\kappa$ as an independent
parameter to be measured from the entanglement entropy rather than using the
theoretical prediction. We note that we cannot measure $\kappa$ from looking
at the correlation length scaling with respect to the bond dimension as the
algorithm at zero temperature already spontaneously breaks translational
symmetry and we lose access to the transfer matrix. Applying the same
methodology with the square CTMRG and assuming the central charge $c=1$ we
found $\eta=0.4839\pm 0.0131$. Furthermore, assuming $\eta=1/2$ we obtain
$c=1.0334\pm 0.0280$. The critical exponent and central charge obtained from
the square CTMRG are less accurate than the ones obtained with the
$C_{3}$-symmetric CTMRG with larger error bars. This is due to noise present
in the fitting of the order parameter and entanglement entropy seen in Fig.
16, and in that instance the $C_{3}$-symmetric CTMRG outperforms the square
CTMRG.

Figure 16: Simulations done at $h=0$ and $T=0$ with the square and honeycomb
CTMRG algorithms using $\chi\in[200,400]$. Top panel : the logarithm of the
order parameter versus the logarithm of the bond dimension. Bottom panel :
entanglement entropy versus the logarithm of the bond dimension. Results
obtained from the honeycomb CTMRG are more accurate than the one obtained from
the square CTMRG. The black lines indicate the slope predicted by Eq. 23 using
$c=1$ while the yellow lines are the linear fits from which we extract the
critical exponents.

We also note that the exponential divergence of the correlation length limits
the range of numerically reachable temperatures. Indeed, since the numerically
accessible correlation length is also bound by the bond dimension with
$\xi_{\chi}=\chi^{\kappa}$, at low enough temperature the (physical)
correlation length will be larger than what can be described with finite bond
dimension $\chi$, $i.e.$ $\xi>\xi_{\chi}$ and the system will freeze in the
sense that (i) the correlation length saturates to a finite value and (ii)
correspondingly, the environment converges to a temperature-independent fixed
point. Due to the exponential divergence of the correlation length, in order
to overcome this problem and to effectively describe the lower temperatures we
would need to increase $\chi$ exponentially fast as well.

###  IV.2 Constrained model

We now turn to the constrained model by setting $K^{-1}=0$. As already
mentioned, there are two phases separated by a KT transition characterized by
$\eta_{KT}=4/9$. We thus use the critical value of $\eta$ as the criterion to
locate the transition and we found $H_{KT}=0.305\pm 0.006$. The exponent is
computed with the scaling relations discussed in the previous section and the
results are shown in Fig. 17 (upper panel). For $H\leq 0.2$ we used bond
dimension $\chi\in[200,400]$ while for $H\geq 0.25$, closer to the transition,
we used $\chi\in[300,520]$. We have plotted the transfer matrix results
obtained from finite-size scaling assuming power law correction (Ref. 83,
Table II) as well. At low reduced fields, our results agree reasonably well
with the transfer matrix results while at larger reduced field, the two
methods start to give different results. It is worth noting that although the
values of the exponent obtained from transfer matrix above $H>0.3$ are
attributed in Ref. 83 to strong cross-overs into the ordered phase (see also
[89, 105]), using the $\eta_{KT}=4/9$ criteria to determine the location of
the KT transition would lead to a larger value than the present result.

Indeed, around the transition, significant corrections to scaling are expected
[83] and the critical field $H_{KT}=0.266\pm 0.01$ in Refs. 82, 83, 84 is only
obtained by considering _both_ finite size logarithmic and power-law
corrections. It is unclear how to perform similar corrections in the framework
of finite-bond dimension scaling. We note that, close to $H\sim 0.3$, our
results show a slight dependence on the range of bond dimension over which we
fit: when considering higher bond dimensions, the computation of $\eta$ gives
a lower exponent. We illustrate the dependency in Fig. 17 (bottom panel) for
three different values of the reduced field where in the $x$-axis $\chi$ is
used as an abuse of notation referring to the interval
$\mathcal{I}_{\chi}=[\chi-250,\chi]$. As expected, the difference between
$\eta_{A}$ and $\eta_{B}$ decreases upon increasing the bond dimension. At
$H=0$ and $H=0.25$ the CTMRG results are well converged. In contrast, at
$H=0.29$ we observe some dependency and it becomes harder to conclude. From
the behavior of $\eta$ in the lower-panel of Fig. 17, we draw the conclusion
that the obtained critical field should be taken as an upper bound for the
transition.

Figure 17: Upper panel : Simulations done at $K^{-1}=0$. At $H=0$, the results
are in good agreement with the exact results $\eta=1/2$. For $H\leq 0.2,$ we
used bond dimension $\chi\in[200,400]$, while for $H>0.2$ we used
$\chi\in[300,520]$. The transfer matrix results are from Ref. 83. Lower panel
: simulations done at $K^{-1}=0$. The $\circ$ and $\triangle$ symbols denote
$\eta_{A}$ and $\eta_{B}$, respectively. In the $\chi$ goes to infinity limit
one recovers the exact result.

###  IV.3 High field

At high field, the critical temperature is large and both transfer matrices
and Monte Carlo methods also have no problem identifying the location of the
Potts transition or its universality class[78, 79, 84]. As an additional
benchmark of our approach, we verify the nature of the transition at $h=3$: if
the transition belongs to the three-state Potts universality class, both
$\psi^{9}$ and $\xi^{-6/5}$ should behave linearly and intersect the $x$-axis
at the critical temperature. This is indeed what we observe in Fig. 18. More
precisely, the intersection between the linear fit of $\psi^{9}$ and the
$x$-axis gives a critical temperature $T_{c}=1.3440$ while the intersection
between the linear fit of $\xi^{-6/5}$ and the $x$-axis gives a critical
temperature $T_{c}=1.3437$. The difference between the two temperatures is of
the order of $10^{-4}$, in agreement with a unique transition belonging to the
three-state Potts universality class. The correlation length has been
extrapolated using bond dimension ranging from $\chi=100$ to $\chi=200$. We
have used $\chi=200$ and $\chi=250$ for the order parameter and we can see in
Fig. 18 that it has converged with respect to the bond dimension.

Figure 18: Simulations were done at $h=3$. The order parameter has been
computed at fixed bond dimension. The critical temperatures obtained from the
order parameter and correlation length respectively are $T_{c}=1.3440$, and
$T_{c}=1.3437$.

###  IV.4 Low field

We now discuss in detail the results obtained at low temperatures. We observe
a strong bond dimension dependency for the order parameter and entanglement
entropy, even away from the critical line. Yet, as previously discussed, such
dependency can be used to determine the critical temperatures and critical
exponents. In contrast to the previous section, we do not work at fixed field
to check the nature of the universality class, but we are looking for the
location of the critical reduced field knowing that the transition is in the
three-state Potts universality class. We thus choose to locate the transition
by looking at the regime where the order parameter decays with respect to
$\chi$ as a power law and perform a self-consistent check by assuming the
central charge $c=4/5$ and then measuring the spin-spin decay critical
exponent $\eta$.

We show the results in Fig. 20 for the lowest considered temperature
$K^{-1}=0.25$ where we used bond dimension up to $\chi=400$. The algebraic
decay of the order parameter can be seen in only a narrow interval
$H_{c}\in[0.40,0.42]$. By considering higher bond dimensions, the interval’s
width would diminish. But we are limited to finite bond dimension and thus
consider $H_{c}=0.41\pm 0.01$. Assuming $c=4/5$ and considering only the
largest bond dimensions between $\chi=200$ and $\chi=400$ we find
$\eta_{A}=0.270\pm 0.003$ and $\eta_{B}=0.266\pm 0.003$ in good agreement with
the three-state Potts universality class. By recovering the right critical
exponent we confirm that we are not in a cross-over regime influenced by the
proximity of the KT point and that we indeed have located the critical reduced
field.

Figure 19: Order parameter and entanglement entropy at temperature
$K^{-1}=0.4$ for different bond dimensions. Figure 20: log-log plot of the
order parameter with respect to the bond dimension at $K^{-1}=0.25$. We can
identify the power law regime to be in between $H=0.40$ and $H=0.42$. We thus
estimate $H_{c}=0.41\pm 0.01$. At $H_{c}=0.41$, by fitting the order parameter
for $\chi\in[200,400]$ one finds $\eta_{A}=0.270\pm 0.003$ and
$\eta_{B}=0.266\pm 0.003$ in agreement with the three-state Potts universality
class.

###  IV.5 Phase Diagram

The full phase diagram is given in Fig. 21, where we have plotted the critical
temperature obtained from the $C_{3}$ symmetric honeycomb CTMRG as well as the
transfer matrix results obtained in Ref. 84. At finite temperature, our
results agree reasonably well with the data from the transfer matrix methods,
but with significantly smaller error bars. Furthermore, we are able to reach
lower temperatures while keeping reasonably small errorbars, such that the
transition line we found shows the first signs of a curvature, which could not
be observed previously. In that regard, our study further supports the
renormalization group predictions which suggest that the Potts critical line
meets the zero temperature KT point with a leading exponent one-half and a
linear correction [70, 84]:

$K_{c}^{-1}\propto(H_{c}-H_{KT})^{1/2}+o(H_{c}-H_{KT}).$ (28)

We have also fitted the new transition line up to the second order in that
expansion and show the results in Fig. 21. We find that the fit of the
transition line is perfectly compatible with both critical reduced fields
$H_{KT}=0.266\pm 0.01$ and $H_{KT}=0.305\pm 0.006$ and is not enough to
distinguish between the two.

As previously mentioned, due to the exponential divergence of the correlation
length upon approaching the $K^{-1}=0$ critical line, there is a limitation
for the temperatures we could reach, and we were not able to study the phase
diagram for non zero temperature smaller than $K^{-1}<0.25$.

Figure 21: Phase diagram of the triangular antiferromagnetic Ising model in a
field shown as $K^{-1}$ with respect to the reduced field $H$. Our results are
in good agreement with previous transfer matrix study. Blue and red dashed
lines indicate the fits done up to second order on the present transition
considering respectively $H_{KT}=0.266\pm 0.01$ and $H_{KT}=0.305\pm 0.006$.

##  V Discussion and Summary

We have introduced a variant of the CTMRG algorithm that contracts infinite
honeycomb tensor networks. It is especially powerful for problems that have a
natural $C_{3}$-symmetry. In the case of the antiferromagnetic triangular
Ising model, this method is more accurate than the square CTMRG and when
considering an extra field outperforms the more standard transfer matrix
methods. By exploring the triangular Ising antiferromagnet in a longitudinal
field, we confirmed that the tensor network construction based on ground-state
rules [53] also allows to reach low temperatures in the honeycomb lattice
CTMRG. Our results mostly confirm and support the predictions from the
literature [70, 83, 84]: at finite temperature, as expected, we found the
nature of the transition to be three-state Potts by measuring a unique
critical temperature from the correlation length and order parameter but also
by using bond dimension scaling of the entanglement entropy and order
parameter. Overall, we confirm Qian et al. results by recovering essentially
the same critical temperatures but with significantly smaller error bars, and
we extend the transition line to lower temperatures previously not available.
As the temperature is lowered to zero, we found evidence of the reduced
critical field converging to the Kosterlitz-Thouless critical reduced field of
the constrained model in favor of the scenario proposed by Nienhuis et al.. At
$T=0$, we locate the KT transition at a slightly higher reduced field than
previously reported.

Several improvements of this CTMRG for the honeycomb lattice might be worth
investigating. The choice of isometry could be different: the most commonly
used isometry in the square lattice CTMRG is the one first introduced by
Corboz for contracting iPEPS [91]. It would be interesting to formulate the
equivalent isometry for CTMRG on honeycomb and to see whether it leads to a
better convergence rate for some problems. A more challenging improvement
would be to design a more general multisite version of the algorithm which can
take arbitrary unit cell; as mentioned, this could be useful to have access to
the transfer matrix also in presence of symmetry breaking. Furthermore, this
algorithm can also be used to contract two-dimensional quantum systems iPEPS
wave-functions on the honeycomb lattice [61] and can easily be combined with
simple or full updates. It could also be implemented with automatic
differentiation, in a similar spirit as Ref. 66. In the framework of PESS or
iPESO, the algorithm could also be used to investigate 2D quantum systems on
the kagome or triangular lattices [106, 107, 108, 109].

Acknowledgments. We thank Loic Herviou, Olivier Gauthé, Mithilesh Nayak,
Andrej Gendiar and Bram Vanhecke for useful discussions. This work has been
supported by the Swiss National Science Foundation grant No. 212082. The
calculations have been performed using the facilities of the Scientific IT and
Application Support Center of EPFL.

## References

  * Nishino [1995] T. Nishino, Density matrix renormalization group method for 2d classical models, J. Phys. Soc. Jpn. 64, 3598 (1995).
  * Nishino and Okunishi [1996] T. Nishino and K. Okunishi, Corner transfer matrix renormalization group method, J. Phys. Soc. Jpn. 65, 891 (1996), https://doi.org/10.1143/JPSJ.65.891 .
  * Nishino and Okunishi [1997] T. Nishino and K. Okunishi, Corner transfer matrix algorithm for classical renormalization group, J. Phys. Soc. Jpn. 66, 3040 (1997).
  * Levin and Nave [2007] M. Levin and C. P. Nave, Tensor renormalization group approach to two-dimensional classical lattice models, Phys. Rev. Lett. 99, 120601 (2007).
  * Orús and Vidal [2009] R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009).
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part ii, Phys. Rev. 60, 263 (1941).
  * Onsager [1944] L. Onsager, Crystal statistics. I. A two-dimensional model with an order-disorder transition, Phys. Rev. 65, 117 (1944).
  * Kaufman [1949] B. Kaufman, Crystal statistics. ii. partition function evaluated by spinor analysis, Phys. Rev. 76, 1232 (1949).
  * Syôzi [1951] I. Syôzi, Statistics of kagomé Lattice, Prog. Theor. Phys. 6, 306 (1951).
  * Wannier [1950] G. H. Wannier, Antiferromagnetism. the triangular Ising net, Phys. Rev. 79, 357 (1950).
  * Houtappel [1950] R. Houtappel, Order-disorder in hexagonal lattices, Physica 16, 425 (1950).
  * Fisher [1966] M. E. Fisher, On the dimer solution of planar Ising models, J. Math. Phys. 7, 1776 (1966).
  * Evenbly and Vidal [2015] G. Evenbly and G. Vidal, Tensor network renormalization, Phys. Rev. Lett. 115, 180405 (2015).
  * Vidal [2007] G. Vidal, Classical simulation of infinite-size quantum lattice systems in one spatial dimension, Phys. Rev. Lett. 98, 070201 (2007).
  * Jordan  _et al._ [2008] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Phys. Rev. Lett. 101, 250602 (2008).
  * Orús and Vidal [2008] R. Orús and G. Vidal, Infinite time-evolving block decimation algorithm beyond unitary evolution, Phys. Rev. B 78, 155117 (2008).
  * Zauner-Stauber  _et al._ [2018] V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Variational optimization algorithms for uniform matrix product states, Phys. Rev. B 97, 045145 (2018).
  * Fishman  _et al._ [2018] M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for contracting infinite two-dimensional tensor networks, Phys. Rev. B 98, 235148 (2018).
  * Nietner  _et al._ [2020] A. Nietner, B. Vanhecke, F. Verstraete, J. Eisert, and L. Vanderstraeten, Efficient variational contraction of two-dimensional tensor networks with a non-trivial unit cell, Quantum 4, 328 (2020).
  * Xie  _et al._ [2009] Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, and T. Xiang, Second renormalization of tensor-network states, Phys. Rev. Lett. 103, 160601 (2009).
  * Zhao  _et al._ [2010] H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai, and T. Xiang, Renormalization of tensor-network states, Phys. Rev. B 81, 174411 (2010).
  * Okunishi  _et al._ [2022] K. Okunishi, T. Nishino, and H. Ueda, Developments in the tensor network — from statistical mechanics to quantum entanglement, Journal of the Physical Society of Japan 91, 062001 (2022).
  * Baxter [1968] R. J. Baxter, Dimers on a rectangular lattice, Journal of Mathematical Physics 9, 650 (1968), https://doi.org/10.1063/1.1664623 .
  * Baxter [1978] R. J. Baxter, Variational approximations for square lattice models in statistical mechanics, J Stat Phys 19, 461 (1978).
  * Tsang [1979] S. K. Tsang, Square lattice variational approximations applied to the Ising model, J Stat Phys 20, 95 (1979).
  * Baxter [1981] R. J. Baxter, Corner transfer matrices, Physica A: Statistical Mechanics and its Applications 106, 18 (1981).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * White [1993] S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions (2004), arXiv:0407.066 [cond-mat.str-el] .
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * Jiang  _et al._ [2008] H. C. Jiang, Z. Y. Weng, and T. Xiang, Accurate determination of tensor network state of quantum lattice models in two dimensions, Phys. Rev. Lett. 101, 090603 (2008).
  * Phien  _et al._ [2015] H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing, Phys. Rev. B 92, 035142 (2015).
  * Liao  _et al._ [2019] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Differentiable programming tensor networks, Phys. Rev. X 9, 031041 (2019).
  * Hasik  _et al._ [2021] J. Hasik, D. Poilblanc, and F. Becca, Investigation of the Néel phase of the frustrated Heisenberg antiferromagnet by differentiable symmetric tensor networks, SciPost Phys. 10, 012 (2021).
  * Corboz  _et al._ [2010] P. Corboz, R. Orús, B. Bauer, and G. Vidal, Simulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair states, Phys. Rev. B 81, 165104 (2010).
  * Corboz  _et al._ [2014a] P. Corboz, T. M. Rice, and M. Troyer, Competing states in the $t$-$j$ model: Uniform $d$-wave state versus stripe state, Phys. Rev. Lett. 113, 046402 (2014a).
  * Corboz  _et al._ [2011] P. Corboz, S. R. White, G. Vidal, and M. Troyer, Stripes in the two-dimensional $t$-$j$ model with infinite projected entangled-pair states, Phys. Rev. B 84, 041108(R) (2011).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Stripe order in the underdoped region of the two-dimensional hubbard model, Science 358, 1155 (2017).
  * Murg  _et al._ [2009] V. Murg, F. Verstraete, and J. I. Cirac, Exploring frustrated spin systems using projected entangled pair states, Phys. Rev. B 79, 195119 (2009).
  * Corboz and Mila [2013] P. Corboz and F. Mila, Tensor network study of the shastry-sutherland model in zero magnetic field, Phys. Rev. B 87, 115144 (2013).
  * Hasik [2019] J. Hasik,  _Towards next-generation methods to optimize two-dimensional tensor networks: Algorithmic differentiation and applications to quantum magnets_ , Ph.D. thesis, SISSA, Trieste (2019).
  * Schmoll  _et al._ [2023] P. Schmoll, A. Kshetrimayum, J. Naumann, J. Eisert, and Y. Iqbal, Tensor network study of the spin-$\frac{1}{2}$ heisenberg antiferromagnet on the shuriken lattice, Phys. Rev. B 107, 064406 (2023).
  * Vanderstraeten  _et al._ [2018] L. Vanderstraeten, B. Vanhecke, and F. Verstraete, Residual entropies for three-dimensional frustrated spin systems with tensor networks, Phys. Rev. E 98, 042145 (2018).
  * Ueda  _et al._ [2020] H. Ueda, K. Okunishi, K. Harada, R. Krčmár, A. Gendiar, S. Yunoki, and T. Nishino, Finite-$m$ scaling analysis of Berezinskii-Kosterlitz-Thouless phase transitions and entanglement spectrum for the six-state clock model, Phys. Rev. E 101, 062111 (2020).
  * Nyckees  _et al._ [2021] S. Nyckees, J. Colbois, and F. Mila, Identifying the huse-fisher universality class of the three-state chiral potts model, Nuclear Physics B 965, 115365 (2021).
  * Schmoll  _et al._ [2021] P. Schmoll, A. Kshetrimayum, J. Eisert, R. Orús, and M. Rizzi, The classical two-dimensional Heisenberg model revisited: An $SU(2)$-symmetric tensor network study, SciPost Phys. 11, 098 (2021).
  * Nyckees and Mila [2022] S. Nyckees and F. Mila, Tensor network investigation of the hard-square model, Phys. Rev. B 106, 174433 (2022).
  * Vieijra  _et al._ [2021] T. Vieijra, J. Haegeman, F. Verstraete, and L. Vanderstraeten, Direct sampling of projected entangled-pair states, Phys. Rev. B 104, 235141 (2021).
  * Frías-Pérez  _et al._ [2023] M. Frías-Pérez, M. Mariën, D. P. García, M. C. Bañuls, and S. Iblisdir, Collective monte carlo updates through tensor network renormalization (2023), arXiv:2104.13264 [cond-mat.stat-mech] .
  * Li and Yang [2021] H. Li and L.-P. Yang, Tensor network simulation for the frustrated ${J}_{1}\text{$-$}{J}_{2}$ Ising model on the square lattice, Phys. Rev. E 104, 024118 (2021).
  * Wang  _et al._ [2014] C. Wang, S.-M. Qin, and H.-J. Zhou, Topologically invariant tensor renormalization group method for the edwards-anderson spin glasses model, Phys. Rev. B 90, 174201 (2014).
  * Zhu and Katzgraber [2019] Z. Zhu and H. G. Katzgraber, Do tensor renormalization group methods work for frustrated spin systems? (2019), arXiv:1903.07721 [cond-mat.dis-nn] .
  * Vanhecke  _et al._ [2021] B. Vanhecke, J. Colbois, L. Vanderstraeten, F. Verstraete, and F. Mila, Solving frustrated Ising models using tensor networks, Phys. Rev. Res. 3, 013041 (2021).
  * Liu  _et al._ [2021] J.-G. Liu, L. Wang, and P. Zhang, Tropical tensor network for ground states of spin glasses, Phys. Rev. Lett. 126, 090506 (2021).
  * Colbois  _et al._ [2022] J. Colbois, B. Vanhecke, L. Vanderstraeten, A. Smerald, F. Verstraete, and F. Mila, Partial lifting of degeneracy in the ${J}_{1}\text{$-$}{J}_{2}\text{$-$}{J}_{3}$ Ising antiferromagnet on the kagome lattice, Phys. Rev. B 106, 174403 (2022).
  * Song and Zhang [2022] F.-F. Song and G.-M. Zhang, Tensor network approach to the two-dimensional fully frustrated xy model and a chiral ordered phase, Phys. Rev. B 105, 134516 (2022).
  * Song and Zhang [2023] F.-F. Song and G.-M. Zhang, Tensor network approach to the fully frustrated xy model on a kagome lattice with a fractional vortex-antivortex pairing transition, Phys. Rev. B 108, 014424 (2023).
  * Song  _et al._ [2023] F.-F. Song, T.-Y. Lin, and G.-M. Zhang, Unified tensor network theory for frustrated classical spin models in two dimensions (2023), arXiv:2309.05321 [cond-mat.str-el] .
  * Akimenko [2023] S. S. Akimenko, Tensor network construction for lattice gas models: Hard-core and triangular lattice models, Phys. Rev. E 107, 054116 (2023).
  * Mangazeev  _et al._ [2010] V. V. Mangazeev, M. Y. Dudalev, V. V. Bazhanov, and M. T. Batchelor, Scaling and universality in the two-dimensional Ising model with a magnetic field, Phys. Rev. E 81, 060103(R) (2010).
  * Corboz  _et al._ [2012a] P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Spin-orbital quantum liquid on the honeycomb lattice, Phys. Rev. X 2, 041013 (2012a).
  * Corboz  _et al._ [2012b] P. Corboz, K. Penc, F. Mila, and A. M. Läuchli, Simplex solids in su($n$) heisenberg models on the kagome and checkerboard lattices, Phys. Rev. B 86, 041106(R) (2012b).
  * Jahromi  _et al._ [2016] S. S. Jahromi, M. Kargarian, S. F. Masoudi, and A. Langari, Topological spin liquids in the ruby lattice with anisotropic kitaev interactions, Phys. Rev. B 94, 125145 (2016).
  * Jahromi  _et al._ [2018] S. S. Jahromi, R. Orús, M. Kargarian, and A. Langari, Infinite projected entangled-pair state algorithm for ruby and triangle-honeycomb lattices, Phys. Rev. B 97, 115161 (2018).
  * Bauer  _et al._ [2012] B. Bauer, P. Corboz, A. M. Läuchli, L. Messio, K. Penc, M. Troyer, and F. Mila, Three-sublattice order in the su(3) heisenberg model on the square and triangular lattice, Phys. Rev. B 85, 125116 (2012).
  * Lukin and Sotnikov [2023] I. V. Lukin and A. G. Sotnikov, Variational optimization of tensor-network states with the honeycomb-lattice corner transfer matrix, Phys. Rev. B 107, 054424 (2023).
  * Vanderstraeten  _et al._ [2022] L. Vanderstraeten, L. Burgelman, B. Ponsioen, M. Van Damme, B. Vanhecke, P. Corboz, J. Haegeman, and F. Verstraete, Variational methods for contracting projected entangled-pair states, Phys. Rev. B 105, 195140 (2022).
  * Stephenson [2004] J. Stephenson, Ising‐Model Spin Correlations on the Triangular Lattice, Journal of Mathematical Physics 5, 1009 (2004).
  * Nienhuis [1982] B. Nienhuis, Exact critical point and critical exponents of $\mathrm{O}(n)$ models in two dimensions, Phys. Rev. Lett. 49, 1062 (1982).
  * Nienhuis  _et al._ [1984] B. Nienhuis, H. J. Hilhorst, and H. W. J. Blote, Triangular sos models and cubic-crystal shapes, Journal of Physics A: Mathematical and General 17, 3559 (1984).
  * Nienhuis [1987] B. Nienhuis,  _Coulomb gas formulations of two-dimensional phase transitions_ , edited by C. Domb and J. Lebowitz, Phase transitions and critical phenomena, Vol. 11 (Academic Press, London, 1987) p. 1–53.
  * Blote and Hilborst [1982] H. W. J. Blote and H. J. Hilborst, Roughening transitions and the zero-temperature triangular Ising antiferromagnet, Journal of Physics A: Mathematical and General 15, L631 (1982).
  * Henley [2010] C. L. Henley, The “Coulomb phase” in frustrated systems, Annual Review of Condensed Matter Physics 1, 179 (2010).
  * Baxter and Tsang [1980] R. J. Baxter and S. K. Tsang, Entropy of hard hexagons, Journal of Physics A: Mathematical and General 13, 1023 (1980).
  * Metcalf and Yang [1978] B. D. Metcalf and C. P. Yang, Degeneracy of antiferromagnetic Ising lattices at critical magnetic field and zero temperature, Phys. Rev. B 18, 2304 (1978).
  * Stephenson [1970] J. Stephenson, Ising‐Model Spin Correlations on the Triangular Lattice. III. Isotropic Antiferromagnetic Lattice, Journal of Mathematical Physics 11, 413 (1970), https://pubs.aip.org/aip/jmp/article-pdf/11/2/413/7519770/413_1_online.pdf .
  * Jacobsen and Fogedby [1997] J. L. Jacobsen and H. C. Fogedby, Monte carlo study of correlations near the ground state of the triangular antiferromagnetic Ising model, Physica A: Statistical Mechanics and its Applications 246, 563 (1997).
  * Alexander [1975] S. Alexander, Lattice gas transition of He on Grafoil. A continuous transition with cubic terms, Physics Letters A 54, 353 (1975).
  * Kinzel and Schick [1981] W. Kinzel and M. Schick, Phenomenological scaling approach to the triangular Ising antiferromagnet, Phys. Rev. B 23, 3435 (1981).
  * Noh and Kim [1992] J. D. Noh and D. Kim, Phase boundary and universality of the triangular lattice antiferromagnetic Ising model, Int. J. Mod. Phys. B 06, 2913 (1992), publisher: World Scientific Publishing Co.
  * Tamashiro and Salinas [1997] M. N. Tamashiro and S. R. Salinas, Bethe-Peierls approximation for the triangular Ising antiferromagnet in a field, Phys. Rev. B 56, 8241 (1997), publisher: American Physical Society.
  * Blöte  _et al._ [1991] H. W. J. Blöte, M. P. Nightingale, X. N. Wu, and A. Hoogland, Field-induced Kosterlitz-Thouless transition in the zero-temperature triangular Ising antiferromagnet, Phys. Rev. B 43, 8751 (1991), publisher: American Physical Society.
  * Blöte and Nightingale [1993] H. W. J. Blöte and M. P. Nightingale, Antiferromagnetic triangular Ising model: Critical behavior of the ground state, Phys. Rev. B 47, 15046 (1993).
  * Qian  _et al._ [2004] X. Qian, M. Wegewijs, and H. W. J. Blöte, Critical frontier of the triangular Ising antiferromagnet in a field, Phys. Rev. E 69, 036127 (2004).
  * Ràcz [1980] Z. Ràcz, Phase boundary of Ising antiferromagnets near $h={H}_{c}$ and $t=0$: Results from hard-core lattice gas calculations, Phys. Rev. B 21, 4012 (1980).
  * Baxter [1982] R. J. Baxter,  _Exactly solved models in statistical mechanics_ (Academic Press, New York, 1982).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * José  _et al._ [1977] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model, Phys. Rev. B 16, 1217 (1977).
  * Otsuka  _et al._ [2006] H. Otsuka, Y. Okabe, and K. Okunishi, Field-induced Berezinskii-Kosterlitz-Thouless transition and string-density plateau in the anisotropic triangular antiferromagnetic Ising model, Phys. Rev. E 73, 035105(R) (2006).
  * Hwang and Kim [2014] C.-O. Hwang and S.-Y. Kim, Field-induced Kosterlitz–Thouless transition in critical triangular-lattice antiferromagnets, Monte Carlo Methods and Applications 20, 217 (2014).
  * Corboz  _et al._ [2014b] P. Corboz, T. M. Rice, and M. Troyer, Competing states in the $t$-$j$ model: Uniform $d$-wave state versus stripe state, Phys. Rev. Lett. 113, 046402 (2014b).
  * Chatelain and Gendiar [2020] C. Chatelain and A. Gendiar, Absence of logarithmic divergence of the entanglement entropies at the phase transitions of a 2d classical hard rod model, The European Physical Journal B 93, 134 (2020).
  * R. Krčmár and Nishino [2020] A. G. R. Krčmár and T. Nishino, Entanglement-entropy study of phase transitions in six-state clock model, Acta Phys. Pol. A (2020).
  * Daniška and Gendiar [2015] M. Daniška and A. Gendiar, Tensor product variational formulation applied to pentagonal lattice, Journal of Physics A: Mathematical and Theoretical 48, 435002 (2015).
  * Daniška and Gendiar [2016] M. Daniška and A. Gendiar, Analysis of quantum spin models on hyperbolic lattices and bethe lattice, Journal of Physics A: Mathematical and Theoretical 49, 145003 (2016).
  * Gendiar  _et al._ [2012] A. Gendiar, R. Krcmar, S. Andergassen, M. Daniška, and T. Nishino, Weak correlation effects in the Ising model on triangular-tiled hyperbolic lattices, Phys. Rev. E 86, 021105 (2012).
  * Nishino  _et al._ [1996] T. Nishino, K. Okunishi, and M. Kikuchi, Numerical renormalization group at criticality, Physics Letters A 213, 69 (1996).
  * Rams  _et al._ [2018] M. M. Rams, P. Czarnik, and L. Cincio, Precise extrapolation of the correlation function asymptotics in uniform tensor network states with application to the bose-hubbard and xxz models, Phys. Rev. X 8, 041033 (2018).
  * Wannier [1973] G. H. Wannier, Antiferromagnetism. the triangular Ising net, Phys. Rev. B 7, 5017 (1973).
  * Villain [1977] J. Villain, Spin glass with non-random interactions, Journal of Physics C: Solid State Physics 10, 1717 (1977).
  * Tanaka  _et al._ [1992] K. Tanaka, T. Horiguchi, and T. Morita, Critical indices for the two-dimensional Ising model with nearest-neighbor and next-nearest-neighbor interactions, Physics Letters A 165, 266 (1992).
  * Pollmann  _et al._ [2009] F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, Theory of finite-entanglement scaling at one-dimensional quantum critical points, Phys. Rev. Lett. 102, 255701 (2009).
  * Tagliacozzo  _et al._ [2008] L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. Latorre, Scaling of entanglement support for matrix product states, Phys. Rev. B 78, 024410 (2008).
  * Schmoll  _et al._ [2019] P. Schmoll, A. Haller, M. Rizzi, and R. Orús, Quantum criticality on a chiral ladder: An su(2) infinite density matrix renormalization group study, Phys. Rev. B 99, 205121 (2019).
  * Hwang and Kim [2010] C.-O. Hwang and S.-Y. Kim, Yang–lee zeros of triangular Ising antiferromagnets, Physica A: Statistical Mechanics and its Applications 389, 5650 (2010).
  * Xie  _et al._ [2014] Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, Tensor renormalization of quantum many-body systems using projected entangled simplex states, Phys. Rev. X 4, 011025 (2014).
  * Schmoll  _et al._ [2022] P. Schmoll, C. Balz, B. Lake, J. Eisert, and A. Kshetrimayum, Finite temperature tensor network algorithm for frustrated two-dimensional quantum materials (2022), arXiv:2211.00121 [cond-mat.str-el] .
  * Jahromi  _et al._ [2020] S. S. Jahromi, R. Orús, D. Poilblanc, and F. Mila, Spin-$\frac{1}{2}$ kagome Heisenberg antiferromagnet with strong breathing anisotropy, SciPost Phys. 9, 092 (2020).
  * Liao  _et al._ [2021] Y. D. Liao, H. Li, Z. Yan, H.-T. Wei, W. Li, Y. Qi, and Z. Y. Meng, Phase diagram of the quantum Ising model on a triangular lattice under external field, Phys. Rev. B 103, 104416 (2021).