# Finite-size and finite bond dimension effects of tensor network
renormalization

Atsushi Ueda<EMAIL_ADDRESS>Institute for Solid State Physics,
University of Tokyo, Kashiwa 277-8581, Japan  Masaki Oshikawa  Institute for
Solid State Physics, University of Tokyo, Kashiwa 277-8581, Japan  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Trans-scale Quantum
Science Institute, University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan

###### Abstract

We propose a general procedure for extracting the running coupling constants
of the underlying field theory of a given classical statistical model on a
two-dimensional lattice, combining tensor network renormalization (TNR) and
the finite-size scaling theory of conformal field theory. By tracking the
coupling constants at each scale, we are able to visualize the renormalization
group (RG) flow and demonstrate it with the classical Ising and 3-state Potts
models. Furthermore, utilizing the new methodology, we reveal the limitations
due to finite bond dimension $D$ on TNR applied to critical systems. We find
that a finite correlation length is imposed by the finite bond dimension in
TNR, and it can be attributed to an emergent relevant perturbation that
respects the symmetries of the system. The correlation length shows the same
power-law dependence on $D$ as the “finite entanglement scaling” of the matrix
product states.

##  I Introduction

The universality of critical phenomena is one of the most intriguing and
important concepts in statistical physics. The renormalization group (RG),
proposed by Wilson [1, 2, 3], provides a conceptual framework to comprehend
and characterize this universality.

In the RG framework, a universality class of critical phenomena is governed by
an RG fixed point in the “theory space”. Theory space is the abstract space of
all possible models or theories that can describe a physical system. Each
point in this space represents a unique combination of the parameters of the
theory, or more concretely, the corresponding Hamiltonian or action. In the
context of the RG, we explore this “theory space” by starting from a specific
point in the theory space and applying the RG transformations. These
transformations effectively move us through the theory space, changing the
values of the parameters as we coarse-grain the system. Importantly, models
within the same universality class converge to an identical position through
the RG transformations. This allows for diverse critical phenomena to be
comprehended in terms of perturbations to the fixed-point Hamiltonian and
their respective scaling behavior.

This theoretical approach has directly facilitated the development of concrete
schemes for the calculation of critical exponents. A prime example is the
$\epsilon$-expansion for the $\phi^{4}$ theory in $4-\epsilon$ dimensions [4,
5]. While the practical utility of such a scheme for calculating critical
exponents may appear to be limited, it is imperative to underscore that the RG
framework establishes the conceptual foundation for understanding the
universality of critical phenomena.

In particular, the fixed point displays conformal invariance in two
dimensions, thereby simplifying the associated theory which is described by a
conformal field theory (CFT). The effective Hamiltonian near the RG fixed
point Hamiltonian, denoted $\hat{H}_{\text{CFT}}$, can be expressed as
follows:

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{\text{CFT}}+\sum_{j}g_{j}\int_{0}^{L}dx\;\hat{\Phi}_{j}(x),$
(1)

In this expression, $\hat{\Phi}_{j}(x)$ represents a scaling operator with a
scaling dimension $x_{j}$, and $g_{j}$ denotes the corresponding coupling
constant. In two dimensions, the running coupling constants $g_{j}$ are
renormalized as $g_{j}\propto L^{2-x_{j}}$ as functions of a scale
$l=\ln(L/a)$. In general, there are only a few RG-relevant coupling constants
with $x_{j}<2$, which increase as the scale $l$ increases.

There also exists RG-irrelevant coupling constants with $x_{j}>2$ that
decrease as the scale $l$ increases. Despite their occasional significance,
the principal characteristics of critical phenomena can be outlined primarily
by considering the limited number of RG-relevant coupling constants.
Differential equations, termed RG equations, frequently describe the evolution
of these running coupling constants as functions of the scale $l$. Field
theory methods frequently serve as the basis for deriving these RG equations.

However, it is worth noting that the exact determination of RG equations may
not always be feasible when the corresponding field theory is not solvable.
Moreover, the application of RG to lattice models has generally been
challenging for quantitative calculations. While offering an intuitive
understanding of RG, the “block spin transformation” method falls short as a
practical computational method for lattice models. Overall, early RG schemes
for lattice models saw limited success, the notable exception being Wilson’s
numerical renormalization group for impurity problem [3].

Subsequently, density matrix renormalization group(DMRG) [6] emerged as a
highly effective numerical algorithm for one-dimensional quantum many-body
systems. Despite its name, DMRG is typically employed as a numerical algorithm
with less emphasis on RG flows in the “theory space.”

More recently, the development of tensor network renormalization(TNR) [7, 8,
9, 10, 11, 12] opened a way to implement numerical schemes for a wide range of
lattice models in a manner more faithful to the original concept of RG.
Notably, it is possible to obtain a fixed-point tensor after multiple
iterations of TNR steps. This fixed-point tensor encapsulates critical
information about the infrared (IR) fixed point, including conformal data.

Regarding RG flow, there have been numerous previous studies[13, 8, 10, 11,
14]. Yet, a generic and quantitative framework for calculating RG flows
remains elusive, primarily due to challenges in maintaining the correlation
between the changes in numerically obtained tensor networks and the RG flow
within the ’theory space.’

In this paper, we first propose an efficient and quantitative scheme to
extract the RG flow numerically from TNR, discussed in Sec. III. This involves
comparing the finite-size spectrum of the transfer matrix with CFT. Concrete
examples, such as the numerical results of the Ising and 3-state Potts models,
are employed to validate the theoretical predictions. Our method also provides
an efficient and accurate estimation of the critical point, extending the
’Level Spectroscopy’ technique previously developed for Berezinskii-
Kosterlitz-Thouless (BKT) transitions [15, 16].

Leveraging this methodology, we uncover the effects of finite bond dimension
$D$ on TNR at criticality in Section IV. The finite-bond approximation of
tensors constrains the effective correlation length, preventing the attainment
of a ’true fixed point tensor’ corresponding to a nontrivial RG fixed point
through repeated TNR procedures. While this phenomenon was reported in earlier
studies [7, 8, 9, 10, 11, 12], it has been often overlooked. Our numerical
results suggest that the finite bond-dimension effects can be attributed to an
emergent relevant perturbation that respects the symmetry of the lattice
model. Furthermore, we demonstrate that the finite correlation length that is
imposed by the finite bond dimension scales in the same way as in matrix
product states (MPS).

We note that some of the methods and observations discussed in this paper were
previously reported in our earlier publication [16], where they were applied
to the classical XY model. The goal of the present paper is to illustrate the
more widespread applicability of this approach and deliver a more
comprehensive analysis of the finite bond-dimension effects.

Sample codes necessary to reproduce the figures presented in this paper, along
with introductory reviews on TRG and TNR, are accessible via Jupyter notebooks
at the following GitHub repository: https://github.com/dartsushi/Loop-
TNR_RGflow.

##  II Review on tensor network renormalization and conformal field theory

Figure 1: A schematic picture of the tensor network renormalization. The
effective local Boltzmann weight at $n$-th RG step $\mathcal{T}^{(n)}$ is
decomposed into the two three-leg tensors and recombined as
$\mathcal{T}^{(n+1)}$. The effective system size enlarges by $\sqrt{2}$ each
RG step. The typical bond dimension and the maximum number of RG steps
employed in this paper are $D\leq 40$, and RG steps $\leq 30$, respectively.

###  II.1 Tensor network renormalization

The tensor network is a numerical technique used to represent the partition
function of statistical models. The partition functions of two-dimensional
statistical models with a system size of $L$ can be expressed through the
contraction of $L^{2}$ tensors. Each tensor represents a local Boltzmann
weight, and its dimensions correspond to physical degrees of freedom. For
instance, the local tensor of the Ising model on the square lattice is a four-
leg tensor
$\mathcal{T}^{(1)}_{ijkl}=e^{\beta(s_{i}s_{j}+s_{j}s_{k}+s_{k}s_{l}+s_{l}s_{i})}$.
The tensor network representation often provides an efficient method for
simulating complex systems.
However, the exact contraction of $L^{2}$ tensors is generally impracticable
for larger system sizes due to the constraints imposed by the high-dimensional
Hilbert space. TNR aims to circumvent this issue by utilizing the principles
of renormalization group theory. During each step of the RG process,
$\mathcal{T}^{(n)}$ is coarse-grained to $\mathcal{T}^{(n+1)}$ via a series of
decompositions and recombinations, as illustrated in Fig. 1. Starting from the
local tensor $\mathcal{T}^{(1)}$, we can simulate a system size of
$L=\sqrt{2}^{n}$ after $n$ RG steps. The process of coarse-graining in TNR
involves numerical truncation, reducing the number of degrees of freedom while
preserving essential physics. Consequently, TNR facilitates efficient
numerical simulation of complex systems.

###  II.2 Critical phenomena under conformal invariance

The following sections mainly discuss the Ising and 3-state Potts models on
the square lattice. The energy (classical Hamiltonian) of the Ising and
3-state Potts models are

$\displaystyle\mathcal{E}_{\rm Ising}$ $\displaystyle=-\sum_{\langle
i,j\rangle}s_{i}s_{j}-h\sum_{i}s_{i},$ (2) $\displaystyle\mathcal{E}_{\rm
Potts}$ $\displaystyle=-\sum_{\langle i,j\rangle}\delta_{s_{i},s_{j}},$ (3)

where $s_{i}=\pm 1$(Ising) and $s_{i}=0,1,2$(3-state Potts). The first terms
and $h$ represent the nearest-neighbor interactions and the magnetic field.
Employing the temperature $T$, the Boltzmann weight is defined as
$e^{-\mathcal{E}/T}$, where we set the Boltzmann constant to unity. Our
primary focus in the main text is the Ising model, while a detailed discussion
of the 3-state Potts model is provided in the appendix. The Ising model
reaches its critical point at $(T,h)=(T_{c},0)$, where
$T_{c}=2/\ln{(1+\sqrt{2})}$. At this criticality, physical quantities like the
spin-spin correlation function are governed by the Ising CFT, which comprises
three primary operators: the identity operator $I$, magnetic operator
$\sigma$, and energy operator $\epsilon$.

In the context of the lattice model, a shift from the critical temperature and
the application of a magnetic field correspond to the perturbative insertion
of $\epsilon$ and $\sigma$ into the effective Hamiltonian. As a result,
$\sigma$ is odd in the $\mathbb{Z}_{2}$ spin-flip, while $I$ and $\epsilon$
are even. Given the operator structure of the CFT, certain quantities are
consequently fixed.

The two-point correlation function is defined as

$\langle\Phi_{i}(r_{i})\Phi_{j}(r_{j})\rangle=\frac{\delta_{i,j}}{|r_{i}-r_{j}|^{2x_{i}}},$

where $\Phi_{i}$ represents a primary operator, and $x_{I}=0$,
$x_{\sigma}=\frac{1}{8}$, and $x_{\epsilon}=1$ are the scaling dimensions. In
a similar vein, the three-point correlation function adopts a universal form,
represented as

$\langle\Phi_{i}(r_{i})\Phi_{j}(r_{j})\Phi_{k}(r_{k})\rangle=\frac{C_{ijk}}{|r_{i}-r_{j}|^{\Delta_{ij}^{k}}|r_{j}-r_{k}|^{\Delta_{jk}^{i}}|r_{k}-r_{i}|^{\Delta_{ki}^{j}}},$

where $C_{ijk}$ is an operator product expansion(OPE) coefficient, and
$\Delta_{ij}^{k}=x_{i}+x_{j}-x_{k}$. The non-zero OPE coefficients are given
as follows:

$\displaystyle C_{III}$
$\displaystyle=C_{I\sigma\sigma}=C_{I\epsilon\epsilon}=1,$ (4) $\displaystyle
C_{\sigma\sigma\epsilon}$ $\displaystyle=\frac{1}{2}.$ (5)

The permutation of the indices does not change the OPE coefficients111The
combination of the indices in non-zero $C_{ijk}$ preserves $\mathbb{Z}_{2}$
symmetry..(For further details of CFT, we suggest readers see Ref. [18].) The
collection of information on the scaling dimension $x_{i}$ and $C_{ijk}$,
referred to as the CFT data, is crucial to understanding critical phenomena.
As such, determining the CFT from a numerical standpoint is of paramount
importance.

##  III Computation of field-theory data and RG flow from TNR

###  III.1 Scaling dimensions

For simplicity, let us consider a classical statistical model on the square
lattice with nearest-neighbor interactions only. Then the local Boltzmann
weight can be represented by a tensor $\mathcal{T}$ with four open indices,
and the partition function is given by contraction of a tensor network which
consists of the tensor $\mathcal{T}$.

More specifically, the partition function $Z(L_{x},L_{y})$ for the system of
the size $L_{x}\times L_{y}$ is given by the contraction of the network of
$L_{x}\times L_{y}$ identical tensors $\mathcal{T}$.

Under a single step of TNR, the length scale represented by a single tensor is
renormalized by $\sqrt{2}$. After $N$ steps, the renormalized tensor becomes
$\mathcal{T}(L)$, which represents the length scale $L=\sqrt{2}^{N}$. As
$\mathcal{T}(L)$ is equivalent to the $L\times L$ contracted tensor network up
to truncation errors, contractions of the horizontal and vertical legs yield
the partition function $Z(L,L)$ in periodic boundary condition (PBC).
Similarly, contracting only the legs in the $x$-direction gives the $L$-stacks
of the transfer matrix in $y$-direction. Since one can regard the transfer
matrix as the imaginary-time evolution operator of corresponding one-
dimensional quantum systems, its eigenvalues $\lambda_{n}(L)$ are related to
the energy levels $E_{n}(L)$ of the quantum system as
$\lambda_{n}(L)=\exp\left(-LE_{n}(L)\right)$. For convenience, we define the
rescaled energy levels $x_{n}(L)$ by $E_{n}(L)-E_{0}(L)=2\pi x_{n}(L)/L$222In
the classical systems, the characteristic velocity $v$, playing the role of
the speed of light, is unity because the system invariant under the exchange
of the $x$ and $y$ axes. to obtain

$\displaystyle\frac{\lambda_{n}(L)}{\lambda_{0}(L)}=\exp(-2\pi x_{n}(L)).$ (6)

Exactly at the criticality, this rescaled energy level $x_{n}(L)$ coincides
with the scaling dimension $x_{n}$ of the in the thermodynamic limit
$(L\rightarrow\infty)$ [20, 21, 22].

If the system is off-critical and without a spontaneous symmetry breaking, the
rescaled energy level of the “first excitation” is asymptotically proportional
to the system size as $\Delta L/(2\pi)$, where $\Delta$ is the inverse
correlation length (excitation gap) in the thermodynamic limit.

Summarizing these observations, naively speaking, we can judge whether the
system is critical or not by looking at the asymptotic behavior of the
rescaled energy levels $x_{n}(L)$. If they grow linearly in $L$, the system is
off-critical. If they approach to constants, the system is critical, and the
scaling dimensions of the operators can be read off from the asymptotic values
of $x_{n}(L)$ in the thermodynamic limit. While this can be a useful guide,
there are corrections from RG-irrelevant perturbations, and more importantly,
due to the limitation of a finite bond dimension, as we will discuss later.

###  III.2 Operator product expansion coefficients

Operator product expansion is another fundamental concept in field theory and
statistical mechanics [23, 24]. Since OPE coefficients determine the structure
of the field theory, their computation is quite important. Numerical
computation of OPE coefficients [25, 26] has not been so straightforward
compared to that of scaling dimensions. Here, we present a simpler way to
compute them, which is applicable to TRG [7], HOTRG [27], and Loop-TNR [10].

As explained in the previous section, the renormalized tensor
$\mathcal{T}^{(N)}$ contracted in $x$-direction is a transfer matrix in the
$y$-direction. While the eigenvalues of the transfer matrix correspond to the
energy or scaling dimension of the primary operators, the eigenvectors thereof
are the wavefunctions of the corresponding “primary states”
$|\psi_{n}(L)\rangle$. This is graphically represented below.

Note that the tensor has been rotated for ease of viewing. We do not change
the contracted index. Likewise, we can compute the wavefunctions of the system
size $2L$ as depicted below.

$|\psi_{n}(L)\rangle$ and $|\psi_{n}(2L)\rangle$ are one-leg and two-leg
tensors, respectively. Thus, one can calculate the overlaps
$\ket{\psi_{\alpha}(2L)}$ and
$\ket{\psi_{\beta}(L)}\otimes\ket{\psi_{\gamma}(L)}$ by contracting the
indices.

$\psi_{\alpha}$ | $\psi_{\beta},\psi_{\gamma}$ | $C_{\alpha\beta\gamma}$ | $2^{2\Delta_{\beta}+2\Delta_{\gamma}-\Delta_{\alpha}}{A_{\alpha\beta\gamma}}/{A_{III}}$
---|---|---|---
$I$ | $\sigma,\sigma$ | 1 | 0.8938
$\sigma$ | $\sigma,I$ | 1 | 0.9473
$I$ | $\epsilon,\epsilon$ | 1 | 0.9966
$\epsilon$ | $\epsilon,I$ | 1 | 0.9968
$\epsilon$ | $\sigma,\sigma$ | 0.5 | 0.5007
$\sigma$ | $\sigma,\epsilon$ | 0.5 | 0.2705
Table 1: The numerically obtained OPE coefficients of the Ising CFT from TRG.
The bond dimension and the system size are $D=56$ and $L=16\sqrt{2}$(9 RG
steps), respectively.

In CFT, the overlap
$\langle{\psi_{\alpha}(2L)}|{\psi_{\beta}(L)}{\psi_{\gamma}(L)}\rangle$ is
proportional to the “pants diagram” of path integrals [28, 29, 30], and the
OPE coefficient and the overlap are related as

$\displaystyle\frac{A_{\alpha\beta\gamma}}{A_{III}}=2^{\Delta_{\alpha}-2\Delta_{\beta}-2\Delta_{\gamma}}C_{\alpha\beta\gamma},$
(7)

where $A_{\alpha\beta\gamma}$, $\Delta,$ and $I$ are
$\langle{\psi_{\alpha}(2L)}|{\psi_{\beta}(L)}{\psi_{\gamma}(L)}\rangle$, the
scaling dimension, and the identity operator, respectively. In most cases, the
identity operator corresponds to the ground state. We benchmark our method by
the critical Ising model. Table. 1 shows the numerically obtained OPE
coefficients by TRG [7] at $L=16\sqrt{2}$ and $D=56$. Naturally, there are
finite-size corrections to Eq. (7). Since Eq. (7) is exact in the
thermodynamic limit, using a very large system size $L$ might appear
desirable. However, as we will discuss later in Sec. IV, corrections due to
the finite bond-dimension effect appear for system sizes larger than a
correlation length $\xi(D)$333This effect is even stronger and non-trivial for
TRG. As reported in Ref. [30], the finite-size effects are significant for
$C_{\sigma\sigma\epsilon}$ and $C_{\epsilon\epsilon I}$. Nevertheless, even
with the moderate size $L=16\sqrt{2}$, the obtained values
$C_{I\epsilon\epsilon}=0.9966$ and $C_{\epsilon\sigma\sigma}=0.5007$ are
rather close to exact CFT results. While we tested our method by the simplest
algorithm, Levin and Nave’s TRG, the method for calculating OPE is
straightforwardly applicable to other TRG and TNR algorithms, such as HOTRG
[27].

###  III.3 Level Spectroscopy

As we have mentioned earlier, the rescaled energy levels $x_{n}(L)$ in Eq. (6)
would be independent of the scale $L$ and give the scaling dimensions of the
corresponding operators, if the system were exactly described by a CFT.
However, the rescaled energy levels of a lattice model generally depend on the
system size $L$, as the effective Hamiltonian of the system contains
perturbations to the CFT as in Eq. (1).

The rescaled energy levels in a finite-size perturbed CFT are given as [20,
21]

$\displaystyle x_{n}(L)=x_{n}+2\pi\sum_{j}C_{nnj}g_{j}(L),$ (8)

where $g_{j}(L)$ scales as $\propto L^{2-x_{j}}$. Comparing Eq. (6) from TNR
and Eq. (8) from the conformal perturbation theory, we can obtain the running
coupling constants $g_{j}(L)$ at each scale from the finite-size effect
$\delta x_{n}(L)=x_{n}(L)-x_{n}$.

An immediate application of this observation is an accurate determination of
the critical point. While such a framework is dubbed “level spectroscopy” was
developed for BKT transition, which is notoriously difficult for standard
finite-size scaling analysis, first for quantum spin systems in one dimension
[15] and recently extended for classical statistical systems in two dimensions
using TNR [16], the basic idea is also applicable to more conventional
critical phenomena such as in the Ising model.

The RG fixed point for the two-dimensional Ising model has two relevant
operators, the energy density $\epsilon$ and the magnetization density
$\sigma$. The coupling constant $g_{\epsilon}$ for $\epsilon$ is proportional
to the deviation of the temperature from the critical point, and also scaled
$\sim L$ in the small coupling limit $g_{\epsilon}\ll 1$ because
$x_{\epsilon}=1$. Thus

$\displaystyle g_{\epsilon}(L)\sim\alpha(T-T_{c})L,$ (9)

when $g_{\epsilon}(L)\ll 1$. Likewise, the coupling $g_{\sigma}$ is
proportional to the magnetic field $h$ and scaled $\sim L^{15/8}$ because
$x_{\sigma}=1/8$.

Although the Ising critical phenomena are mostly described by the two relevant
coupling constants $g_{\epsilon}$ and $g_{\sigma}$, more accurate description
can be obtained by including irrelevant perturbations. Including the leading
irrelevant operators, namely the irrelevant operators with the smallest
scaling dimension permitted by the symmetries, we obtain

$\displaystyle H=H^{*}_{\rm Ising}+\int_{0}^{L}dx$
$\displaystyle[g_{\sigma}\sigma(x)+g_{\epsilon}\epsilon(x)$
$\displaystyle+g_{T^{2}}T_{\rm cyl}^{2}(x)+g_{\bar{T}^{2}}\bar{T}_{\rm
cyl}^{2}(x)],$ (10)

where $T_{\text{cyl}}$ and $\bar{T}_{\text{cyl}}$ are the holomorphic and
anti-holomorphic parts of stress tensor on a cylinder [21]. The holomorphic
part $T_{\text{cyl}}$ of the stress tensor on a cylinder is related to that on
the infinite plane $T_{zz}(z)$ via the conformal mapping $z=e^{2\pi w/L}$,
where $w=\tau+ix$ and $0\leq x<L$. More explicitly, $T_{zz}(z)$ transforms as

$\displaystyle T_{\rm cyl}(w)$
$\displaystyle=\left(\frac{2\pi}{L}\right)^{2}\left(z^{2}T_{zz}(z)-\frac{c}{24}\right).$
(11)

This leads to

$\displaystyle T_{\text{cyl}}(x)$
$\displaystyle=\frac{2\pi}{L}\left(\sum_{n=-\infty}^{\infty}L_{n}e^{2\pi
ix/L}-\frac{c}{24}\right),$ (12)

where $c$ is the central charge characterizing the CFT, and $L_{n}$’s are
generators of the Virasoro algebra defined by

$\displaystyle T_{zz}(z)=\sum_{n=-\infty}^{\infty}\frac{L_{n}}{z^{n+2}},$ (13)

in terms of the holomorphic part $T_{zz}$ of the energy-momentum tensor on the
infinite plane. Inserting the above $T_{\text{cyl}}$ and integrating over
$0\leq x<L$ with an appropriate regularization, the $g_{T^{2}}$-term of the
perturbation is given as [32]

$\displaystyle\int dx\;T^{2}_{\rm
cyl}(x)=L_{0}^{2}-\frac{c+2}{12}L_{0}+2\sum_{n=1}^{\infty}L_{-n}L_{n}+\frac{c(22+5c)}{2880}$

Only the first and second terms affect the energy levels, and the
contributions to $x_{\sigma}(L)$ and $x_{\epsilon}(L)$ are calculated to be
$-\frac{7}{768}g_{T^{2}}$ and $\frac{7}{48}g_{T^{2}}$ respectively. The
computation of the contributions from $\bar{T}^{2}$ is exactly the same, and
we denote their sum as $g$. These operators are the leading irrelevant
operators for the Ising model on the square lattice. Although they break the
continuous rotation symmetry (which corresponds to the Lorentz invariance in
the Minkowski space-time), they are allowed on the square lattice, which is
invariant only under the discrete C4 rotation. The calculation of
$x_{\sigma}(L)$ and $x_{\epsilon}(L)$ is straightforward, and they are shown
in Table. 2 444As $T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ are not
primary operators, we need to pay special attention. The details are discussed
in the Appendix..

While the exact critical point is known for the Ising model on the square
lattice, let us demonstrate the determination of the critical point from the
TNR spectrum without using prior knowledge of the critical point (but
utilizing the CFT data, assuming that we identify the universality class).
Since we are interested in the critical point at zero magnetic fields, we can
set $g_{\sigma}\propto h=0$. The simplest way to determine the critical point
is to look at the lowest rescaled energy level $x_{\sigma}(L)$ in the lowest
order of the relevant coupling constant $g_{\epsilon}$, ignoring the
irrelevant perturbation $g$. Within this approximation, the shift $\delta
x_{\sigma}(L)=x_{\sigma}(L)-x_{\sigma}$ vanishes at the critical point
$T=T_{c}$ where $g_{\epsilon}=0$. Away from the critical point, $\delta
x_{\sigma}(L)$ is non-zero and grows proportionally to $L$ because
$g_{\epsilon}(L)$ scales as $L$. Because of this, we can identify the critical
point with the temperature where $\delta x_{\sigma}(L)=0$ is observed in the
TNR spectrum. However, this estimate suffers from the corrections due to the
leading irrelevant perturbations $T_{\text{cyl}}^{2}$ and
$\bar{T}_{\text{cyl}}^{2}$. Since they have scaling dimension $4$, the
corresponding coupling constant is renormalized as $g\propto L^{-2}$. This
leads to an error of $O(L^{-2})$ in the naive estimate of the critical point
using $\delta x_{\sigma}(L)=0$.

Figure 2: Example of estimating the transition temperature using Loop-TNR. We
set $T^{-}=2.66$ and $T^{+}=2.68$ as an initial estimate. The level-crossing
temperature $T^{*}(L)$ is linearly fitted to extrapolate the transition
temperature. The insert shows how we compute $T^{*}(L)$ for various system
sizes.

We can improve the accuracy by removing the effects of the leading irrelevant
perturbation $g$. This can be done by combining the shifts of the rescaled
energy levels $\delta x_{\sigma}(L)$ and $\delta x_{\epsilon}(L)$ as

$\displaystyle\delta x_{\text{cmb}}\equiv$ $\displaystyle\delta
x_{\sigma}(L)+\frac{1}{16}\delta x_{\epsilon}(L)$ $\displaystyle=\pi
g_{\epsilon}+(\alpha^{\sigma}_{\sigma}+\frac{1}{16}\alpha_{\epsilon}^{\sigma})g_{\sigma}^{2}+(\alpha^{\epsilon}_{\sigma}+\frac{1}{16}\alpha^{\epsilon}_{\epsilon})g_{\epsilon}^{2}.$
(14)

Note that the first-order correction in the irrelevant coupling $g$ is
canceled out. Now we can identify the critical point by finding the
temperature for which $\delta x_{\text{cmb}}\propto g_{\epsilon}(L)=0$. Having
eliminated the effects of the leading irrelevant perturbation
$T_{\text{cyl}}^{2},\bar{T}_{\text{cyl}}^{2}$, the dominant error is now
caused by the next-leading irrelevant operator with scaling dimension $6$ and
thus should be scaled as $L^{-4}$.

In practice, the determination of the critical point can be efficiently
implemented as follows. First, we pick up one temperature from each phase:
$T^{+}>T_{c}$ and $T^{-}<T_{c}$, and calculate the combined shift $\delta
x_{\text{cmb}}$ at these temperatures. The phase of the system can be
confirmed by observing the growth of $\delta x_{\text{cmb}}$ as the system
size increases because it increases/decreases if the system is in the high-
temperature/low-temperature phase (if the initial choice of the temperature
turns out to be wrong, change the temperature and restart the process). Next,
linear interpolations of the combined shift between the two temperatures
$T^{\pm}$ are made, and the crossing of the lines for system sizes $L$ and
$\sqrt{2}L$ is found, as shown in the insert of Fig. 2. We denote the
temperature where the two lines cross as $T^{*}(L)$. Because of the second-
order contribution $O({g_{\epsilon}}^{2})$ in Eq. (14), the crossing
temperature $T^{*}(L)$ obtained by the linear interpolation deviates from the
true critical point $T_{c}$ as $T^{*}(L)-T_{c}\propto g_{\epsilon}\propto L$,
when $g_{\epsilon}\ll 1$555It is proportional to $L^{2-x_{\text{thermal}}}$,
where $x_{\text{thermal}}$ is the scaling dimension of the thermal operator..
The critical point $T_{c}$ is estimated by fitting $T^{*}(L)$ by a linear
function of $L$ as $T^{*}(L)\sim T_{c}+\text{const.}L$. While the
“extrapolation” to $L=0$ used here might look unusual, this procedure is done
to remove the effect of the nonlinearity due to $O({g_{\epsilon}}^{2})$ in Eq.
(14), and the condition $\delta x_{\text{cmb}}=0$ itself is accurate for
$T_{c}$ up to the error of $O(L^{-4})$ due to the next-leading irrelevant
perturbations. An example of the estimate of $T_{c}$ with the above procedure
with the choice of the temperatures $T^{+}=2.68$ and $T^{-}=2.66$ and with
system sizes $16\leq L<64$ is depicted in Fig. 2. The final estimate of the
critical point is $T_{c}^{\text{est}}=2.269177$. Remarkably, even with the
choice of two temperatures differ by $10^{-2}$ and the relatively low bond-
dimension $D=20$, the estimated critical point is quite accurate:
$T_{c}^{\text{est}}-T_{c}=-8.11\times 10^{-6}$. This is thanks to the
suppression of the error to $O(L^{-4})$ by eliminating the contributions from
the leading irrelevant operators. Once the critical point is estimated with
good accuracy with this procedure, the accuracy can be further improved by
choosing $T^{\pm}$ closer to the estimated critical temperature and then
applying the same procedure.

model | operator | Rescaled energy level
---|---|---
Ising model |  |
$x_{\sigma}(L)$ | $\frac{1}{8}+\alpha_{\sigma}^{\sigma}g_{\sigma}^{2}+\pi g_{\epsilon}+\alpha^{\epsilon}_{\sigma}g_{\epsilon}^{2}-\frac{7}{768}\pi g$
$x_{\epsilon}(L)$ | $1+\alpha^{\sigma}_{\epsilon}g_{\sigma}^{2}+\alpha^{\epsilon}_{\epsilon}g_{\epsilon}^{2}+\frac{7}{48}\pi g$
Table 2: The finite-size scaling dimension of the Ising model. $\alpha$ is a
constant determined from the second-order perturbation. Since $g_{T^{2}}$ and
$g_{\bar{T}^{2}}$ decay in the same manner, we write them as $g$.

###  III.4 Renormalization Group flow

The comparison between the TNR spectrum (6) and the conformal perturbation
theory (8) can also be used to extract running coupling constants and their
scale dependence, enabling a visualization of the RG flow. While this was
shown for the BKT transition in the XY model [16], here let us demonstrate the
method for the Ising model. This will also be useful to investigate the finite
bond-dimension effects in detail, as we will discuss in Sec. IV. The
extraction of running coupling constants in the Ising model is again based on
the shifts of the rescaled energy levels in Table 2, and it is useful to
consider the combined shift (14) also for this purpose. Given $g_{\sigma}$ and
$g_{\epsilon}$ are small in the vicinity of criticality, we neglect
$g_{\epsilon}^{2}$ for $h=0$ and redefine two relevant coupling constants as
$g_{t}=\pi g_{\epsilon}$ and
$g_{h}=\sqrt{(\alpha^{\sigma}_{\sigma}+\frac{1}{16}\alpha_{\epsilon}^{\sigma})}g_{\sigma}$
for convenience. In this way, the combined shift Eq. (14) simply gives $g_{t}$
when $h=0$ and ${g_{h}}^{2}$ when $T=T_{c}$, in the lowest order of
$g_{t},g_{h}$. Using these relations, we can read off the relevant coupling
constants $g_{t}$ or $g_{h}$ from the TNR data, as shown in Fig. 3($b$).

Figure 3: (Left panel) The system size dependence of $\delta
x_{\text{cmb}}=\delta x_{\sigma}+\delta x_{\epsilon}/16$ for $h=\pm
10^{-5}$(purple and green), $T=1.0001T_{c}$(red) and $T=0.9999T_{c}$(blue).
The purple and green dots are on top of each other, and “$+$” denotes the data
with a negative sign. After removing the $L^{-2}$ irrelevant perturbations,
the next leading $L^{-4}$ perturbation shown with a blue dotted line appears.
The data was obtained via Loop-TNR with a bond dimension of $D=24$, which was
deemed sufficient for the finitely-correlated systems being considered. (Right
panel) The resulting renormalization group flow. Only data after six steps are
exhibited, where the $L^{-4}$ perturbations disappear.

As we have discussed in the previous subsection, the effects of the leading
irrelevant perturbations $T_{\text{cyl}}^{2},\bar{T}_{\text{cyl}}^{2}$ with
scaling dimension $4$ are eliminated in the combined shift (14), and thus the
finite-size correction is now of $O(L^{-4})$, due to the next-leading
irrelevant operators with scaling dimension $6$. This $O(L^{-4})$ scaling is
indeed observed in Fig. 3 near the critical point for small system size $L$
when relevant perturbations are still negligible. Since it is safe to say that
these contributions disappear after five RG steps, we can conclude that the
origin of $g_{t}$ and $g_{h}$ are purely from $\epsilon$ and $\sigma$ after
six steps.

The right panel illustrates the scale-dependence of the coupling constants
$g_{t}$ and $g_{h}$. It is nothing but the RG flow of the Ising critical
point, and we conclude that we succeed in calculating the RG flow of the
celebrated Ising fixed point.

There is one thing to note on the left panel of Fig. 3. While the combined
shift (14), which is an estimator for $|g_{h}|^{2}$, scales as $L^{3.75}$ at
$L<10^{3}$, it starts to flatten and scales as $L$ at $L>10^{3}$. This
behavior has a rather simple origin. Since the magnetic perturbation is
relevant, the system has a finite correlation length or equivalently, a non-
zero gap $\Delta$. This implies that the rescaled energy levels are
proportional to $L$ for sufficiently large system size $L\gg\Delta^{-1}$, as
mentioned in Sec. III.1. As a consequence, the shift (14) also grows
proportionally to $L$. In this regime, the conformal perturbation theory
breaks down (higher-order contributions are important), and we no longer
identify the shift (14) with $|g_{h}|^{2}$. This should be distinguished from
the $L$-linear behavior of the combined shift (14) observed for $L>10$ with
$h=0$ and $T\neq T_{c}$, which corresponds to the renormalization of
$g_{t}\propto L$ because of $x_{\epsilon}=1$. The $L$-linear behavior due to
the gap is observed in the non-perturbative regime $\delta
x_{\epsilon,\sigma}\gg x_{\epsilon,\sigma}$, whereas the $L$-linear behavior
due to the scaling is observed in the perturbative regime $\delta
x_{\epsilon,\sigma}\ll x_{\epsilon,\sigma}$.

##  IV Finite bond-dimension effects

Let us examine the impacts of a finite bond-dimension $D$ on TNR from the
perspective of our method. In any computation that employs tensor networks, it
is necessary to restrict the bond dimension to a finite value $D$ due to the
increasing storage requirements and computational costs associated with larger
bond dimensions. The finiteness of the bond dimension inevitably leads to a
loss of information in each step of renormalization after a certain number of
iterations. Although TNR can nominally handle arbitrary large systems, and the
TNR-type calculations are often used to study extremely large systems, we have
to be careful about the limitations due to the finite bond dimension.

The limitation of the finite bond dimension $D$ on the MPS is characterized by
the finite (maximum) correlation length $\xi(D)$ of the MPS [35, 36, 37]. The
correlation length of MPS is known to obey the scaling law

$\displaystyle\xi(D)\sim$ $\displaystyle D^{\kappa},$ (15)
$\displaystyle\kappa=$ $\displaystyle\frac{6}{c(1+\sqrt{\frac{12}{c}})}.$ (16)

While the TNR-type calculation of two-dimensional statistical systems appears
rather different from the MPS applied to one-dimensional quantum systems, the
emergence of the finite correlation length $\xi(D)$ obeying the similar
scaling law (15) was reported in Ref. [38] for a HOTRG calculation of the
critical Ising model in two dimensions. The exponent $\kappa$ for the Ising
model was estimated to be approximately $2$, which is close to the MPS
exponent (16) $\kappa=2.03425\ldots$ for the Ising CFT with central charge
$c=1/2$. A similar emergence of the finite correlation length $\xi(D)$ was
also reported in our TNR finite-size scaling study of the two-dimensional XY
model [16], with the MPS exponent (16) for $c=1$.

In the following, using our TNR finite-size scaling methodology, we will
demonstrate that the emergence of the finite correlation length due to the
finite bond dimension in TNR can be attributed to an emergent relevant
perturbation (Sec. IV.1). Furthermore, we present evidences for the scaling
(15) with the MPS exponent (16) in TNR of Ising and 3-state Potts models (Sec.
IV.2).

###  IV.1 Emergent relevant perturbation

Figure 4: Shift $|\delta x_{\sigma}(L)|$ for the Ising model at $T=T_{c}$,
$h=0$ computed by Loop-TNR with $D=32$. There is little finite-$D$ effect for
small system sizes $L<256$. The emergent perturbations of $\epsilon$ and
$\sigma$ appear at $L\sim 256$ and $L\sim 10^{4}$, scaling as $L$ and
$L^{15/4}$. The induced gap by finite-$D$ goes towards constant at $L>10^{5}$
as denoted with the purple dotted line.

If a finite correlation length emerges in the TNR, it would be natural to
identify the renormalized tensor with a Hamiltonian for the system away from
the critical point, that is, an RG fixed-point (CFT) Hamiltonian perturbed
with relevant operators

$\displaystyle
H_{\text{FB}}(D)=H^{*}_{\text{CFT}}+\sum_{i}\int_{0}^{L}dxg_{i}(D,L)\Phi_{i}(x,D),$
(17)

where $H_{FB}$ is the effective Hamiltonian of the finite-$D$ system and
$\Phi_{i}(x,D)$ are the scaling operators representing the perturbations. In
this view, we expect relevant perturbations to emerge in order to mimic the
finite correlation length imposed by the finite bond dimension.

To demonstrate the emergence of the relevant perturbation, we investigate the
system-size dependence of the shift in the rescaled energy levels $\delta
x_{\sigma}$. In Fig. 4, we show the absolute value of the shift $|\delta
x_{\sigma}|$ as a function of the system size $L$ used in calculating the
transfer matrix spectrum in TNR exactly at the critical point $h=0,T=T_{c}$.
The conformal perturbation theory in Eq. (8) implies that the shift
$x_{\sigma}$ contains contributions from the irrelevant perturbations. Since
the leading irrelevant operators at the critical points are
$T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ with scaling dimension
$4$, we expect $\delta x_{\sigma}(L)$ decays as $L^{-2}$. (This is to be
contrasted with Eq. (14) and Fig. 3, in which the contributions from
$T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ are eliminated.) The
expected $L^{-2}$ behavior in the shift $\delta x_{\sigma}(L)$ is indeed
observed for small system sizes $L<256$. For larger system sizes, however,
$|\delta x_{\sigma}(L)|$ starts to increase, deviating from the conformal
perturbation theory scaling $L^{-2}$. We identify the finite bond-dimension
$D$ effects as the origin of this deviation. More remarkably, we can observe a
clear scaling behavior of the deviation. That is, the shift $|\delta
x_{\sigma}(L)|$ scales with the system sizes as $L$ and $L^{15/4}$ for
$256<L<10^{4}$ and $10^{4}<L$, respectively. Compared with the off-critical
cases in Fig. 3, we realize that these scalings are identical to those induced
by the thermal and magnetic perturbations. In other words, the relevant
perturbations emerge in the TNR calculation.

Let us first discuss the $L^{15/4}$ scaling of the shift, observed for
$L>10^{4}$. This can be understood as the effect of an emerging magnetic
perturbation $h$. Although the magnetic perturbation $h$ is forbidden by the
$\mathbb{Z}_{2}$ spin-flip symmetry, the symmetry could be broken by the
limitations in the machine precision. Once the spin-flip symmetry is broken,
the magnetic field $h$, which is a relevant perturbation, is effectively
generated. Even if the effective magnetic field $h$ is extremely small, it
will be enhanced at each RG step and eventually dominates the system at
sufficiently large length scales. This is what we observe for $L>10^{4}$. This
phenomenon should be related to machine precision and not intrinsic to the
algorithm. If we are interested in a $\mathbb{Z}_{2}$ symmetric system, we can
impose the symmetry at each step of TNR in order to avoid this effect.

In contrast, the $L$ scaling observed for $256<L<10^{4}$ is more intrinsic.
The most relevant perturbation allowed under the $\mathbb{Z}_{2}$ symmetry to
the critical Ising fixed point is the thermal operator. Thus, we expect that
the finite bond dimension effect can be mimicked by the thermal perturbation
$\epsilon$ to the fixed-point Hamiltonian $H^{*}_{\text{CFT}}$. If this is the
case, the effective coefficient $g_{\epsilon}$ grows proportionally to $L$ as
the system size $L$ is increased, because the thermal operator $\epsilon$ has
the scaling dimension $1$. According to Eq. (8), this will lead to a
correction proportional to $L$ in the rescaled energy level $\delta
x_{\sigma}(L)$. This is indeed supported by the numerical result shown in Fig.
4.

In general, the finite-$D$ effect in TNR would be described in terms of the
emergence of relevant perturbation(s) to the fixed-point Hamiltonian, which
induces the finite correlation length $\xi(D)$. In addition to the emergence
of the relevant operator $\epsilon$ in the critical Ising model discussed
above, a similar emergence of the relevant operator is observed in the
critical 3-state Potts model, as demonstrated in Appendix A.

###  IV.2 Scaling of the emergent correlation length

Figure 5:  $(a)$ The scaling of the shift $\delta x_{\sigma}$ in TNR of the
Ising model at the critical point, for various bond dimensions
$D=4,\ldots,28$. The vertical axis is scaled as $L^{2}\delta x_{\sigma}$ so
that it is constant when $L\ll\xi(D)$. When $L\ll\xi(D)$, the shift is
dominated by the emergent relevant perturbation $\epsilon$; this is confirmed
by the scaling $L^{2}g_{\epsilon}\propto L^{3}$. The horizontal axis is scaled
as $L/\xi(D)$, where the correlation length $\xi(D)$ is hypothesized as in
Eqs. (15) and (16). The collapse of the data for different bond dimensions is
strong evidence of the hypothesized scaling of the correlation length
$\xi(D)$. The blue dotted line indicates $L/\xi(D)=1$. We set
$\xi(D)=2.0D^{\kappa}$ so that $L/\xi(D)=1$ becomes the crossover scale
between the finite-size scaling regime and the finite-$D$ scaling regime.
$(b)$ Similar scaling analysis of the shift $\delta x_{\epsilon}$ in TNR of
the 3-state Potts model at the critical point, for various bond dimensions
$D=16,\ldots,40$ with $\xi(D)=0.067D^{\kappa}$. The scaled shift
$L^{0.8}\delta x_{\epsilon}$ behaves as a constant in the finite-size scaling
regime $L/\xi(D)<1$, whereas it scales as $L^{3.2}$ in the finite-$D$ scaling
regime $L/\xi(D)>1$, as expected from the CFT analysis (see Appendix A for
details). The data for different bond dimensions collapse again, giving
compelling evidence for the scaling of the correlation length (15) and (16)

Now let us demonstrate that the finite correlation length $\xi(D)$ induced by
the finite bond dimension $D$ in TNR obeys the same scaling (15) and (16) as
in the MPS, as suggested in Refs. [38, 16].

In Fig. 5, we demonstrate the scaling of the correlation length induced by the
finite bond dimension in TNR of the critical Ising and the 3-state Potts
models. In Fig. 5(a), we plot the shift $\delta x_{\sigma}$ in the Ising model
obtained by the TNR of the Ising model at the critical point, which was also
studied in Fig. 4, with the several different bond dimensions $D=4,\ldots,28$.
Here, we rescaled the vertical axis as $L^{2}\delta x_{\sigma}$ so that the
constant behavior is observed for system size smaller than the correlation
length, where the leading irrelevant perturbation (which causes $\delta
x_{\sigma}\propto L^{-2}$) is dominant. The deviation from the constant at
larger system sizes $L$ can be attributed to the emergent relevant
perturbation $\epsilon$ induced by the finite bond dimension $D$, as discussed
in the previous subsection. This is confirmed by the $L^{3}$ scaling ($L^{2}$
times $\delta x_{\sigma}\propto g_{\epsilon}\propto L$). Most importantly, the
horizontal axis is the rescaled system size $L/\xi(D)$ using the hypothesized
correlation length $\xi(D)=aD^{\kappa}$ given by Eqs. (15) and (16). The
collapse of the data for different bond dimensions strongly supports our
hypothesis on the correlation length. Note that we roughly fit the prefactor
$a$ so that the cross-over occurs at $L=\xi(D).$

In order to confirm the finite-$D$ scaling of the correlation length and its
universality, we have also studied the 3-state Potts model at the critical
point. As an example, in Fig. 5(b), we plot the shift of the rescaled energy
level corresponding to the energy operator $\epsilon$ in the 3-state Potts
model. For this shift $\delta x_{\epsilon}$, the contribution from the leading
irrelevant operator is $\sim L^{-4/5}$, and the dominant contribution from the
emergent relevant perturbation $\epsilon$ is expected to be proportional to
${g_{\epsilon}}^{2}\propto L^{12/5}$. (See Appendix A for details). We
rescaled the vertical axis as $L^{0.8}\delta_{\epsilon}$ so that it is
constant in the finite-size scaling regime $L<\xi(D)$. The horizontal axis is
again the rescaled system size $L/\xi(D)$, with the correlation length
$\xi(D)$ defined in Eqs. (15) and (16) with the central charge $c=4/5$ for the
3-state Potts model. The data for different bond dimensions again show a
collapse, providing compelling evidence for our hypothesis on the correlation
length scaling. For $L/\xi(D)>1$, the data fits well the expected behavior
$L^{0.8}\times{g_{\epsilon}}^{2}\propto L^{0.8}\times L^{2.4}=L^{3.2}$.

##  V Conclusion and discussion

In the first part of the paper, we discussed a method for computing the
coupling constants using renormalized tensors based on the finite-size scaling
theory of CFT. By plotting the resulting values at each scale, we were able to
visualize the RG flow, and we confirmed that the theoretical RG flows, as
shown in Fig. 3, are consistent with the Ising and XY models. Our method has
the advantage of being able to extract both ultraviolet and infrared
information, making it a valuable tool for investigating gapped and crossover
systems.

In the second part of the paper, applying the methodology developed in the
first part, we explored the impact of finite bond-dimension $D$ on the RG
flow. The finiteness of the bond-dimension results in a finite correlation
length $\xi(D)$, or equivalently in a non-zero gap in the energy spectrum of
the corresponding one-dimensional quantum system. We find that this gap
formation can be attributed to the emergence of a relevant perturbation
enforced by the finite bond dimension. This is demonstrated by the RG flow of
the emergent relevant coupling.

The finite-size scaling of TNR shows a crossover at $L\sim\xi(D)$, above which
the system is governed by the finite correlation length. The correlation
length $\xi(D)$ induced by the finite bond dimension in TNR shows the same
scaling (15), (16) as the correlation length of MPS. While such scaling in TNR
was suggested earlier in Refs. [38, 16], in this paper, we presented more
convincing evidence.

Although we do not have a mathematical proof for the scaling of $\xi(D)$ in
TNR at this point, it may be natural from the following point of view. Besides
the construction of the transfer matrix by contracting horizontal legs, the
renormalized tensor obtained in TNR can give the corner transfer matrix by
contracting the upper and left legs. The same finite-$D$ scaling (15), (16) as
in MPS was observed in corner transfer matrix renormalization group (CTMRG)
[39, 40, 41]. Moreover, the entanglement spectrum for the half-bipartition of
the system of length $2L$ can be related to a contraction of four renormalized
tensors of linear size $L$ [42], as shown in Fig. 6. These relations are
suggestive of the identical scaling of $\xi(D)$ in MPS, CTMRG, and TNR as we
have observed.

Figure 6:  (Left panel) A schematic picture of the reduced density matrix
$\rho_{A}$ for a bipartition of the system in the path integral picture. The
uncontracted legs correspond to the indices of the reduced density matrix.
(Right panel) Each of the four quadrants of the space-time in the left panel
may be replaced by the renormalized tensor in TNR with appropriate boundary
conditions.

Our study highlights the importance of considering the impact of the finite
bond dimension in the TNR-type approach. In particular, a direct study of the
thermodynamic limit with TNR would be prone to errors due to the finite
correlation length $\xi(D)$ imposed by the finite bond dimension. As a
resolution of this problem, we have demonstrated that accurate data for the
thermodynamic limit can be extracted by finite-size scaling of TNR spectra
obtained for system sizes smaller than $\xi(D)$, combined with conformal field
theory. Even with this limitation, the tractable system size is greatly
increased from $\sim\log{D}$ with exact diagonalization to $\xi(D)\sim
D^{\kappa}$ in TNR.

Note added When this work was almost completed, a closely related work [43]
based on a HOTRG study of the Ising model appeared. It is quite similar in
spirit to this work, combining finite-size scaling and HOTRG. Their estimate
of the correlation length $\xi$ from the transfer matrix eigenvalue, and the
determination of the critical point based on the finite-size scaling of $\xi$,
are essentially equivalent to our analysis of $\delta x_{\sigma}$ discussed in
Sec. III.3. Utilizing the knowledge of Ising CFT, we have further improved the
accuracy by analyzing $\delta x_{\text{cmb}}$ which removes the effects of the
leading irrelevant operators. The error in the estimated transition
temperature $T_{c}^{\text{est}}$ we obtained in Sec. III.3 is about $10^{-5}$,
which is larger than theirs ($10^{-7}\sim 10^{-6}$). However, our estimate is
based on the data at two temperatures $T^{\pm}$ separated by $10^{-2}$ and can
be further improved by taking more data points. The MPS scaling (15) and (16)
of the correlation length in TNR we have discussed is also supported by them.

## Acknowledgements

A. Ueda thanks Tsuyoshi Okubo and Luca Tagliacozzo for stimulating
discussions. This work was supported in part by MEXT/JSPS KAKENHI Grant Nos.
JP17H06462 and JP19H01808, and JST CREST Grant No. JPMJCR19T2, and was
partially done during the program “Tensor Networks: Mathematical Structures
and Novel Algorithms” held in Erwin Schrödinger International Institute for
Mathematics and Physics (ESI) at the University of Vienna. A part of the
computation in this work has been done using the facilities of the
Supercomputer Center, the Institute for Solid State Physics, the University of
Tokyo.

## References

  * Wilson [1971a] K. G. Wilson, Renormalization group and critical phenomena. I. renormalization group and the Kadanoff scaling picture, Physical review B 4, 3174 (1971a).
  * Wilson [1971b] K. G. Wilson, Renormalization group and critical phenomena. II. phase-space cell analysis of critical behavior, Physical Review B 4, 3184 (1971b).
  * Wilson [1975] K. G. Wilson, The renormalization group: Critical phenomena and the Kondo problem, Reviews of modern physics 47, 773 (1975).
  * Fisher [1974] M. E. Fisher, The renormalization group in the theory of critical behavior, Rev. Mod. Phys. 46, 597 (1974).
  * Wilson and Fisher [1972] K. G. Wilson and M. E. Fisher, Critical exponents in 3.99 dimensions, Physical Review Letters 28, 240 (1972).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * 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).
  * Evenbly and Vidal [2015] G. Evenbly and G. Vidal, Tensor network renormalization, Phys. Rev. Lett. 115, 180405 (2015).
  * Evenbly [2017] G. Evenbly, Algorithms for tensor network renormalization, Phys. Rev. B 95, 045117 (2017).
  * Yang  _et al._ [2017] S. Yang, Z.-C. Gu, and X.-G. Wen, Loop optimization for tensor network renormalization, Phys. Rev. Lett. 118, 110504 (2017).
  * Bal  _et al._ [2017] M. Bal, M. Mariën, J. Haegeman, and F. Verstraete, Renormalization group flows of hamiltonians using tensor networks, Phys. Rev. Lett. 118, 250602 (2017).
  * Hauru  _et al._ [2018] M. Hauru, C. Delcamp, and S. Mizera, Renormalization of tensor networks using graph-independent local truncations, Phys. Rev. B 97, 045111 (2018).
  * Zou  _et al._ [2018] Y. Zou, A. Milsted, and G. Vidal, Conformal data and renormalization group flow in critical quantum spin chains using periodic uniform matrix product states, Phys. Rev. Lett. 121, 230402 (2018).
  * Delcamp and Tilloy [2020] C. Delcamp and A. Tilloy, Computing the renormalization group flow of two-dimensional ${\phi}^{4}$ theory with tensor networks, Phys. Rev. Res. 2, 033278 (2020).
  * Nomura and Okamoto [1994] K. Nomura and K. Okamoto, Critical properties of s= 1/2 antiferromagnetic XXZ chain with next-nearest-neighbour interactions, Journal of Physics A: Mathematical and General 27, 5773 (1994).
  * Ueda and Oshikawa [2021] A. Ueda and M. Oshikawa, Resolving the Berezinskii-Kosterlitz-Thouless transition in the two-dimensional XY model with tensor-network-based level spectroscopy, Phys. Rev. B 104, 165132 (2021).
  * Note [1] The combination of the indices in non-zero $C_{ijk}$ preserves $\mathbb{Z}_{2}$ symmetry.
  * Francesco  _et al._ [2012] P. Francesco, P. Mathieu, and D. Sénéchal,  _Conformal field theory_ (Springer Science & Business Media, 2012).
  * Note [2] In the classical systems, the characteristic velocity $v$, playing the role of the speed of light, is unity because the system invariant under the exchange of the $x$ and $y$ axes.
  * Cardy [1984] J. L. Cardy, Conformal invariance and universality in finite-size scaling, Journal of Physics A: Mathematical and General 17, L385 (1984).
  * Cardy [1986] J. L. Cardy, Operator content of two-dimensional conformally invariant theories, Nuclear Physics B 270, 186 (1986).
  * Gu and Wen [2009] Z.-C. Gu and X.-G. Wen, Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order, Phys. Rev. B 80, 155131 (2009).
  * Kadanoff [1969] L. P. Kadanoff, Operator algebra and the determination of critical indices, Phys. Rev. Lett. 23, 1430 (1969).
  * Wilson [1969] K. G. Wilson, Non-lagrangian models of current algebra, Phys. Rev. 179, 1499 (1969).
  * Evenbly and Vidal [2016] G. Evenbly and G. Vidal, Local scale transformations on the lattice with tensor network renormalization, Phys. Rev. Lett. 116, 040401 (2016).
  * Li  _et al._ [2022] G. Li, K. H. Pai, and Z.-C. Gu, Tensor-network renormalization approach to the $q$-state clock model, Phys. Rev. Research 4, 023159 (2022).
  * Xie  _et al._ [2012] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Coarse-graining renormalization by higher-order singular value decomposition, Phys. Rev. B 86, 045139 (2012).
  * Zou and Vidal [2022] Y. Zou and G. Vidal, Multiboundary generalization of thermofield double states and their realization in critical quantum spin chains, Phys. Rev. B 105, 125125 (2022).
  * Zou [2022] Y. Zou, Universal information of critical quantum spin chains from wavefunction overlap, Phys. Rev. B 105, 165420 (2022).
  * Liu  _et al._ [2022] Y. Liu, Y. Zou, and S. Ryu, Operator fusion from wavefunction overlaps: Universal finite-size corrections and application to haagerup model (2022).
  * Note [3] This effect is even stronger and non-trivial for TRG.
  * Poghosyan [2019] A. Poghosyan, Shaping lattice through irrelevant perturbation: Ising model, Journal of High Energy Physics 2019, 1 (2019).
  * Note [4] As $T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ are not primary operators, we need to pay special attention. The details are discussed in the Appendix.
  * Note [5] It is proportional to $L^{2-x_{\text{thermal}}}$, where $x_{\text{thermal}}$ is the scaling dimension of the thermal operator.
  * 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).
  * 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, Physical review letters 102, 255701 (2009).
  * Pirvu  _et al._ [2012] B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, Matrix product states for critical spin chains: Finite-size versus finite-entanglement scaling, Phys. Rev. B 86, 075117 (2012).
  * Ueda  _et al._ [2014] H. Ueda, K. Okunishi, and T. Nishino, Doubling of entanglement spectrum in tensor renormalization group, Phys. Rev. B 89, 075116 (2014).
  * Nishino  _et al._ [1996] T. Nishino, K. Okunishi, and M. Kikuchi, Numerical renormalization group at criticality, Physics Letters A 213, 69 (1996).
  * Ueda  _et al._ [2017] H. Ueda, K. Okunishi, R. Krčmár, A. Gendiar, S. Yunoki, and T. Nishino, Critical behavior of the two-dimensional icosahedron model, Phys. Rev. E 96, 062112 (2017).
  * 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).
  * Calabrese and Lefevre [2008] P. Calabrese and A. Lefevre, Entanglement spectrum in one-dimensional systems, Phys. Rev. A 78, 032329 (2008).
  * Huang  _et al._ [2023] C.-Y. Huang, S.-H. Chan, Y.-J. Kao, and P. Chen, Tensor network based finite-size scaling for two-dimensional ising model, Phys. Rev. B 107, 205123 (2023).
  * Dotsenko [1984] V. Dotsenko, Critical behaviour and associated conformal algebra of the z3 potts model, Nuclear Physics B 235, 54 (1984).
  * Dotsenko and Fateev [1984] V. Dotsenko and V. Fateev, Conformal algebra and multipoint correlation functions in 2D statistical models, Nuclear Physics B 240, 312 (1984).
  * Dotsenko and Fateev [1985a] V. Dotsenko and V. Fateev, Four-point correlation functions and the operator algebra in 2d conformal invariant theories with central charge $c\leq 1$, Nuclear Physics B 251, 691 (1985a).
  * Dotsenko and Fateev [1985b] V. Dotsenko and V. Fateev, Operator algebra of two-dimensional conformal theories with central charge $c\leq 1$, Physics Letters B 154, 291 (1985b).
  * Fuchs and Klemm [1989] J. Fuchs and A. Klemm, The computation of the operator algebra in non-diagonal conformal field theories, Annals of Physics 194, 303 (1989).
  * Esterlis  _et al._ [2016] I. Esterlis, A. L. Fitzpatrick, and D. M. Ramirez, Closure of the operator product expansion in the non-unitary bootstrap, Journal of High Energy Physics 2016, 1 (2016).
  * Note [6] For each iteration, the lattice rotates by 45 degrees, and it corresponds to the conformal transformation $w=e^{\frac{i\pi}{4}}z$ on a complex plane. As the irrelevant perturbations $T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ have a conformal spin $4$ and $-4$, they get an additional factor $(e^{\frac{i\pi}{4}})^{4}=(e^{\frac{i\pi}{4}})^{-4}=-1$ for an odd number of steps. We can see this by plotting the data from even steps (original) and odd steps (tilt) separately.

Symbol | Dimension | Meaning
---|---|---
$I$ | 0 | identity
$\epsilon$ | $\frac{2}{5}$ | thermal op.
$\sigma$ | $\frac{1}{15}$ | spin
$X$ | $\frac{7}{5}$ |
$Y$ | ${3}$ |
$Z$ | $\frac{2}{3}$ |
Table 3: A set of primary operators of the 3-state Potts model.

##  Appendix A Finite-Entanglement scaling of the Three-State Potts model.

Figure 7: The size dependence of the (a)$\delta{x}_{\sigma}$ and
(b)$\delta{x}_{\epsilon}$ at $T=0.999995T_{c}$ and $T=1.000005T_{c}$. The pink
and green dotted lines denote $L^{-0.8}$, (a)$L^{1.2}$, and (b)$L^{2.4}$
fittings respectively. For the low-temperature phase, the sign of
$\delta{x}_{\sigma}$ is negative at $L>100$. The dip on the left panel around
$L\sim 10^{2}$ corresponds to the zero point of Eq. (19). (b) The finite-size
effect to the $x_{\epsilon}$ suffers less from
$T_{\text{cyl}}^{2}+\bar{T}_{\text{cyl}}^{2}$ in amplitude. The scaling of Eq.
(20) is clearly observed.

###  A.1 Model

We can further verify the emergence of relevant perturbations by applying it
to the three-state Potts model. It is a natural extension of the Ising model
to the $\mathbb{Z}_{3}$ symmetry, and the Hamiltonian is

$\displaystyle H=-\sum_{\langle i,j\rangle}\delta_{s_{i},s_{j}},$ (18)

where $s_{i}$ takes 0, 1, and $-1$. It has a phase transition of
$\mathbb{Z}_{3}$ symmetry breaking at $T_{c}=1/\log(1+\sqrt{3})$. The critical
theory of the 3-state Potts model is another type of the minimal model
$\mathcal{M}(6,5)$ with $c=\frac{4}{5}$[18, 44]. A set of primary operators
are shown in Table. 3.

As opposed to the Ising model, there are off-diagonal operators as
$\Phi_{\frac{2}{5},\frac{7}{5}},\ \Phi_{\frac{7}{5},\frac{2}{5}}$ and
$\Phi_{3,0},\ \Phi_{0,3}$(currents).

Let us first examine the RG flow in a gapped system. Similar to the Ising
model, the phase transition is identified by spontaneous symmetry breaking.
The high-temperature phase is a trivial phase, whereas the low-temperature
region is $\mathbb{Z}_{3}$ symmetry breaking phase. Thus, the fixed-point
tensor is a stacking of three states with their $\mathbb{Z}_{3}$ charge $0,\
-1$, and 1.

###  A.2 Construction of the effective Hamiltonian

The RG flow can be seen by investigating the scaling dimensions. For instance,
we can take the spin operator $\sigma=\Phi_{\frac{1}{15},\frac{1}{15}}$ and
plot the value of $\delta{x}_{\sigma}=x_{\sigma}(L)-\frac{2}{15}$. Similarly,
as in the Ising model, there is competition between irrelevant and relevant
operators: $X=\Phi_{\frac{7}{5},\frac{7}{5}}$ and
$\epsilon=\Phi_{\frac{2}{5},\frac{2}{5}}$. The thermal operator separates the
$\mathbb{Z}_{3}$ symmetry-breaking phase from the trivial one. The finite-size
corrections of $X$ and $\epsilon$ to $x_{\sigma}$ are $L^{-0.8}$ and
$L^{1.2}$, respectively. The fusion rules are
$\sigma\times\sigma=1+\epsilon+\sigma+X+Y+Z$, $\epsilon\times\epsilon=1+X$,
and $\epsilon\times\sigma=\sigma+Z$. Hence, $\delta{x}_{\sigma}$ has the
following form:

$\displaystyle\delta{x}_{\sigma}=2\pi c_{\sigma\sigma
X}g_{X}\left(\frac{L}{2\pi}\right)^{-0.8}+2\pi
c_{\sigma\sigma\epsilon}g_{\epsilon}\left(\frac{L}{2\pi}\right)^{1.2}.$ (19)

On the other hand, the perturbation of $\epsilon$ appears as a second-order
term for $\delta x_{\epsilon}$ because the fusion rule says
$\epsilon\times\epsilon=1+X$. Consequently, $\delta x_{\epsilon}$ can be
computed as

$\displaystyle\delta{x}_{\epsilon}=2\pi c_{\epsilon\epsilon
X}g_{X}\left(\frac{L}{2\pi}\right)^{-0.8}+\alpha
g^{2}_{\epsilon}\left(\frac{L}{2\pi}\right)^{2.4},$ (20)

where $\alpha$ is a constant determined from the second-order calculation.

Figure 8: $36\delta x_{\sigma}-\delta x_{\epsilon}$ for the high temperature
phase. “+” is used when the sign is negative. The red dotted line denotes the
$L^{-2}$ fitting while the light green one is just a relevant $L^{1.2}$
contribution from $\epsilon$. Loop-TNR rotates the lattice by $\frac{\pi}{4}$
at each RG step, and the tilted system is plotted with the blue dots.

Figure. 7 shows the computed $\delta x_{\sigma}$ by TNR. As expected, it
exhibits the competition between irrelevant and relevant operators. The sign
of $g_{\epsilon}$ is the opposite between two phases, which is a manifest
indication of the RG flow in the opposite direction due to the thermal
operator. $x_{\sigma}$ has doubly degenerate states with $\mathbb{Z}_{3}$
charge $\pm 1$. In the low-temperature phase, these two states flow to
$x_{\sigma}(L)\rightarrow 0$, and the fixed point tensor becomes three-fold
degenerate. As for the irrelevant perturbation, there seems to be a
discrepancy between $\delta x_{\sigma}$ in Fig. 7 and Eq. (19). The data
points are scattered for small system sizes and not precisely on the fitting
lines. This is due to the leading irrelevant operator we have not considered.
We can identify it as $T_{\text{cyl}}^{2}+\bar{T}_{\text{cyl}}^{2}$ as
followings. Just as we did in the left panel of Fig. 3, the contributions from
$g_{X}$ can be eliminated by combining $\delta x_{\sigma}$ and $\delta
x_{\epsilon}$. The OPE coefficients for the 3-state Potts model are known, and
the ratio of the two OPE coefficients is ${c_{\epsilon\epsilon
X}}/{c_{\sigma\sigma X}}=36$[45, 46, 47, 48, 49]. Thus, the origin of the
”scattering” shall be observed by plotting $36\delta x_{\sigma}-\delta
x_{\epsilon}$.

Figure. 8 displays the result for the high-temperature phase. It is now
obvious that the scattering of Fig. 7 comes from the $L^{-2}$ perturbation
denoted with the red dotted line. Also, it has a conformal spin $s$ because it
flips a sign at each step and $s\equiv 4$ (mod 8) 666For each iteration, the
lattice rotates by 45 degrees, and it corresponds to the conformal
transformation $w=e^{\frac{i\pi}{4}}z$ on a complex plane. As the irrelevant
perturbations $T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$ have a
conformal spin $4$ and $-4$, they get an additional factor
$(e^{\frac{i\pi}{4}})^{4}=(e^{\frac{i\pi}{4}})^{-4}=-1$ for an odd number of
steps. We can see this by plotting the data from even steps (original) and odd
steps (tilt) separately.. As a result, we can conclude the irrelevant operator
has the conformal weights as $(h,\bar{h})=(4,0)$ and $(0,4)$, which are
$T_{\text{cyl}}^{2}$ and $\bar{T}_{\text{cyl}}^{2}$. Finally, the effective
Hamiltonian of the critical 3-state Potts model on the square lattice can be
constructed as

$\displaystyle H=H^{*}_{\rm
Potts}+\int_{0}^{L}dx\left[g_{X}\Phi_{\frac{7}{5},\frac{7}{5}}(x)+g_{T}(T_{\text{cyl}}^{2}+\bar{T}_{\text{cyl}}^{2})\right].$
(21)

###  A.3 Finite-Entanglement scaling

At the critical temperature of the Ising model, the finite-$D$ effect proves
to be a perturbation from the thermal operator. Let us verify it for the
critical 3-state Potts model. Due to the irrelevant perturbations from
$T_{\text{cyl}}^{2}+\bar{T}_{\text{cyl}}^{2}$, the finite-$D$ effects are
clearer for $x_{\epsilon}(L)$ as seen in Fig. 7($b$). This is shown in Fig. 5
of the main text. Here, we demonstrate that $\delta x_{\sigma}(L)$ also shows
the universal behavior with $L/\xi(D)$.

Figure 9: Rescaled $\delta x_{\sigma}$ by $\xi(D)=D^{\kappa}$ at the critical
temperature. The resulting data collapse onto a universal function that is
independent of $L/\xi(D)$. If $L/\xi(D)<1$, the system is in the FSS region,
while if $L/\xi(D)\geq 1$, it is in the FES region. In the FES region, the
scaling of the first-order and second-order perturbations are indicated by a
gray and pink line, respectively. $x_{\sigma}$ is computed as an average value
of the first and second excitation energy.

Figure. 9 shows the rescaled correction to $\delta x_{\sigma}(L)$. For
$L>\xi(D)$, the perturbation grows as $L^{2}$ denoted by a gray line, which
means that the emergent perturbation scales as $L^{1.2}$. Compared with Eq.
(19), it is clear that the emergent perturbation is from the thermal operator.
However, as the system size increases, the second-order perturbation becomes
predominant as shown with a pink line. As $\epsilon$ is the most relevant
operator that is permitted by symmetry, it supports our conjecture stated in
the main text.