Finite correlation length scaling with infinite projected entangled pair states at finite temperature

TL;DR

This paper applies finite correlation length scaling to iPEPS thermal states to accurately determine critical temperatures and exponents of 2D quantum phase transitions at finite temperature.

Contribution

It demonstrates the effectiveness of FCLS in extracting critical properties from iPEPS thermal states at finite temperature, extending previous zero-temperature methods.

Findings

FCLS accurately estimates $T_c$ and critical exponents near phase transitions.

Correlation length remains finite at $T_c$ for accessible bond dimensions.

Results agree with Quantum Monte Carlo except near quantum critical points.

Abstract

We study second order finite temperature phase transitions of the 2D quantum Ising and interacting honeycomb fermions models using infinite projected entangled pair states (iPEPS). We obtain an iPEPS thermal state representation by Variational Tensor Network Renormalization (VTNR). We find that at the critical temperature $T_c$ the iPEPS correlation length is finite for the computationally accessible values of the iPEPS bond dimension $D$. Motivated by this observation we investigate the application of Finite Correlation Length Scaling (FCLS), which has been previously used for iPEPS simulations of quantum critical points at $T=0$, to obtain precise values of $T_c$ and the universal critical exponents. We find that in the vicinity of $T_c$ the behavior of observables follows well the one predicted by FCLS. Using FCLS we obtain $T_c$ and the critical exponents in agreement with Quantum Monte Carlo (QMC) results except for couplings close to the quantum critical points where larger bond dimensions are required.

Equations

$$\rho (T)=\frac{1}{Z(T)}{e}^{-H/T},Z(T)=\text{Tr}{e}^{-H/T}.$$

$$\rho (T)={\text{Tr}}_a|\Psi (T)⟩⟨\Psi (T)|.$$

$$|\Psi(T=\infty)\rangle\propto\bigotimes_{i}\big{(}\sum_{s_{i}=1\dots d}|s_{i}s_{i}\rangle\big{)},$$

$$|\Psi (T)⟩\propto e^{-H/(2T)}|\Psi (T=\infty )⟩.$$

$$H=\sum _{j=1,\dots ,m}{H}_{cl}^j.$$

$$e^{-\tau H}=$$

$${\lambda }_j/{\lambda }_0=e^{-({ϵ}_j+i{\varphi }_j)},$$

$$C_{oo}(R)=⟨o_0{o}_R⟩-⟨o_0⟩⟨o_R⟩=\sum _{j>0}{f}_j^{oo}{e}^{-({ϵ}_j+i{\varphi }_j)R}.$$

$$f_j^{oo}=({\varphi }_0|{\mathcal{T}}_o|{\varphi }_j)({\varphi }_j|{\mathcal{T}}_o|{\varphi }_0),$$

$${\xi }_{\mathcal{T}}=1/{ϵ}_1.$$

$$\delta ={ϵ}_2^{mm}-{ϵ}_1^{mm}.$$

$${ϵ}_1^{mm}={ϵ}_1^e+a{\delta }^b,$$

$$\xi =1/{ϵ}_1^e.$$

$$\xi (t=0,D)={\xi }_D.$$

$$m(t,D)={\xi }_D^{-\beta /\nu }\mathcal{M}(t{\xi }_D^{1/\nu }),$$

$$m(t,D,\chi )={\xi }_D^{-\beta /\nu }\bar{\mathcal{M}}(t{\xi }_D^{1/\nu },{\xi }_D/{\xi }_{\chi }),$$

$$m^{\prime }(t,D)={\xi }_D^{-(\beta -1)/\nu }{\mathcal{M}}^{\prime }(t{\xi }_D^{1/\nu }),$$

$$t\frac{m^{\prime }(t,D)}{m(t,D)}=\bar{\mathcal{P}}(t{\xi }_D^{1/\nu }),$$

$${\xi }_D^{1/\nu }\propto \frac{m^{\prime }(t=0,D)}{m(t=0,D)},$$

$$t\frac{m^{\prime}(t,D)}{m(t,D)}=\mathscr{P}\Big{(}t\frac{m^{\prime}(t=0,D)}{m(t=0,D)}\Big{)},$$

$$H=-\sum _{<i,j>}{\sigma }_z^i{\sigma }_z^j+\Gamma \sum _i{\sigma }_x^i,$$

$$m(t,D){\xi }_D^{\beta /\nu }=\mathcal{M}(t{\xi }_D^{1/\nu }),$$

$$m(t=0,D){\xi }_D^{\beta /\nu }=a,$$

$$\log m(t=0,D)=-\beta /\nu \log {\xi }_D+\log a.$$

$$m(t,D){\xi }_D^{\beta /\nu }=\mathcal{M}(t{\xi }_D^{1/\nu }),$$

$$m(t,D)t^{-\beta }=\bar{\mathcal{M}}(t{\xi }_D^{1/\nu }).$$

$$m(t,D)t^{-\beta}=\bar{\mathscr{P}}\Big{(}t\frac{m^{\prime}(t=0,D)}{m(t=0,D)}\Big{)},$$

$$H=-t\sum_{<i,j>}\Big{(}c_{i}^{\dagger}c_{j}+c_{j}^{\dagger}c_{i}\Big{)}+V\sum_{<i,j>}n_{i}n_{j}.$$

$$m=⟨n_A⟩-⟨n_B⟩,$$

Full text

I Introduction

Tensor networks Verstraete et al. (2008); Schöllwock (2011); Orús (2014); Bridgeman and Chubb (2017); Orus (2018) are representations of weakly entangled states obeying an area law of entanglement Hastings (2007); Eisert et al. (2010); Laflorencie (2016). They are a basis for variational numerical methods for strongly correlated quantum many body systems, enabling simulations of fermionic, bosonic and spin models with the same leading computational complexity Pineda et al. (2010); Corboz et al. (2010a); Corboz and Vidal (2009); Kraus et al. (2010); Barthel et al. (2009); Corboz et al. (2010b). The powerful density matrix renormalization group (DMRG) White (1992, 1993) approximates a state of a system by a 1D tensor network called matrix product state (MPS) Affleck et al. (1987); Fannes et al. (1992); Östlund and Rommer (1995). 2D projected entangled pair states (PEPS) Verstraete and Cirac (2004), called also tensor product states Nishino et al. (2001); Nishio et al. (2004), were initially applied as a variational ansatz for 2D ground states Murg et al. (2007); Jordan et al. (2008); Jiang et al. (2008); Corboz et al. (2010b); Kraus et al. (2010) bringing new insights into paradigmatic models of strongly correlated systems (see e. g. Corboz et al. (2014); Zheng et al. (2017); Xie et al. (2014); Liao et al. (2017); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Corboz and Mila (2014)). Recent years brought new applications of PEPS Kshetrimayum et al. (2017); Vanderstraeten et al. (2018); Kennes (2018); Czarnik et al. (2019); Hubig and Cirac (2019) and further progress in the fields of numerical optimization Lubasch et al. (2014); Phien et al. (2015); Corboz (2016); Vanderstraeten et al. (2016) and contraction Xie et al. (2017); Fishman et al. (2018) of PEPS.

Thermal states of 2D local Hamiltonians obey an area law for mutual information, which is reproduced by projected entangled-pair operators (iPEPO) representing thermal states and infinite projected entangled-pair states (iPEPS) representing purifications of thermal states  Wolf et al. (2008), giving motivation to use iPEPS for thermal states simulations Czarnik et al. (2012); Czarnik and Dziarmaga (2014). Recently new methods for simulation of thermal states, based on iPEPS and iPEPO, were proposed Czarnik and Dziarmaga (2015); Czarnik et al. (2016); Czarnik and Dziarmaga (2018); Kshetrimayum et al. (2019); Czarnik et al. (2019). Recent years brought also developments in the field of the closely related direct contraction methods for 3D tensor networks representing partition functions of 2D quantum models Gu et al. (2008); Li et al. (2011); Xie et al. (2012); Ran et al. (2012, 2013) and MPS/MPO based simulations of thermal states of finite width cylinders Bruognolo et al. (2017); Chen et al. (2018a). Some of those methods were already applied to challenging problems Peng et al. (2017); Czarnik et al. (2017); Ran et al. (2018); Chen et al. (2018b, c).

Among demanding problems in the field of 2D strongly correlated systems are finite temperature critical phenomena and in particular finite temperature second order phase transitions. Some of the methods mentioned above were already applied to investigate 2D critical phenomena based on the assumption that large enough $D$ can be obtained to provide results which are converged in $D$ Peng et al. (2017); Czarnik et al. (2017); Chen et al. (2018b); Czarnik et al. (2019), however, for more challenging cases reaching convergence in $D$ will in general be difficult.

Here we demonstrate that even in the case when convergence in $D$ cannot be obtained, it is possible to take finite $D$ effects systematically into account using a Finite Correlation Length Scaling (FCLS) Tagliacozzo et al. (2008); Pollmann et al. (2009); Corboz et al. (2018); Rader and Läuchli (2018). Furthermore, we show that FCLS can be used to obtain critical data for a finite temperature phase transition, i. e. the critical temperature $T_{c}$ and the universal critical exponents.

FCLS, originally called finite entanglement scaling (FES), was first proposed to investigate 1D quantum critical points by infinite MPS (iMPS) Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012). These critical points violate the area law of entanglement Vidal et al. (2003) and as such cannot be represented by finite $D$ iMPS, which have a finite correlation length $\xi_{D}$. It was shown that in the case of the optimal iMPS finite $D$ ground state approximation the finite $D$ modifies observables of the critical state as if the system was finite with the size proportional to $\xi_{D}$ Tagliacozzo et al. (2008); Pollmann et al. (2009). It was also shown that the scaling of the observables with increasing $\xi_{D}$ can be used to determine the critical exponents and the precise location of the critical point similarly as in standard finite size scaling for Quantum Monte Carlo (QMC) simulations Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012).

A similar idea was applied earlier in corner transfer matrix renormalization group (CTMRG) simulations of 2D critical thermal states of classical models. CTMRG approximately contracts a 2D tensor network representing a partition function of a 2D classical system. The approximation introduces an effective correlation length controlled by a refinement parameter $\chi$ of the method Nishino et al. (1996). A scaling ansatz assuming that this correlation length is proportional to the effective system size was introduced to find the critical properties Nishino et al. (1996).

FCLS was recently applied to iPEPS simulations of 2D Lorentz-invariant quantum critical points, i. e. quantum critical points with a linear dispersion relation of low energy excitations Corboz et al. (2018); Rader and Läuchli (2018). It was shown that in such case the optimal finite $D$ iPEPS approximating a critical ground state has a finite correlation length $\xi_{D}$, and that FCLS can be used to determine the critical coupling and the universal critical exponents.

In this paper we simulate second order finite temperature phase transitions for a 2D quantum Ising model and interacting spinless fermions on a honeycomb lattice using Variational Tensor Network Renormlization (VTNR) Czarnik and Dziarmaga (2015); Czarnik et al. (2016). For thermal states at finite $T$ we can expect that the exact state can be represented with a finite bond dimension $D_{exact}$ Corboz et al. (2018). Here we find that at the critical temperature the obtained thermal states have a finite correlation length $\xi_{D}$ for all bond dimensions used in this work suggesting that we are in a regime where $D<D_{exact}$. This motivates us to investigate the possibility to use FCLS also in these cases. In this paper we present benchmark results demonstrating that indeed FCLS can be applied to determine the critical data.

This paper is organized as follows: In Sec. II we introduce a thermal state’s representation by a purification and in Sec. III we describe how to represent such a purification using an iPEPS. In Sec. IV we introduce VTNR which we use to obtain the purification’s iPEPS representation. In Sec. V we describe the CTMRG method which allows us to efficiently contract a 2D tensor network. In Sec. VI we explain how to determine the correlation length of an iPEPS using CTMRG. In Sec. VII we introduce FCLS and in Sec. VIII we describe how to determine the critical temperature $T_{c}$ for a second order phase transition using FCLS. In Secs. IX and X we present the benchmark results for the application of FCLS to simulations of thermal second order phase transitions in the quantum Ising and interacting honeycomb fermions models, respectively. Finally, we provide our conclusions in Sec. XI.

II Purifications of thermal states

A thermal state of a Hamiltonian $H$ for temperature $T$ is given by its thermal density matrix

[TABLE]

Here we consider lattice models for which the Hilbert space $\mathcal{H}$ is a tensor product of Hilbert spaces of individual lattice sites $\mathcal{H}_{i}$ spanned by states $\{|s_{i}\rangle,s_{i}=1\dots d\}$.

We represent $\rho(T)$ by its purification $|\Psi(T)\rangle$ which is a pure state in an enlarged Hilbert space $\tilde{\mathcal{H}}$ created by introducing ancillary degrees of freedom. $\tilde{\mathcal{H}}$ is a tensor product of enlarged Hilbert spaces of individual sites $\tilde{\mathcal{H}}_{i}$, which are spanned by states $\{|s_{i},a_{i}\rangle,s_{i}=1\dots d,a_{i}=1\dots d\}$ with an index $a_{i}$ numbering the ancillary degrees of freedom. To obtain $\rho(T)$ from $|\psi(T)\rangle$ one needs to trace out the ancillary degrees of freedom

[TABLE]

For $T=\infty$ we have

[TABLE]

and for finite $T$ $|\Psi(T)\rangle$ is obtained by an action of $e^{-H/(2T)}$ on the physical degrees of freedom of ${|\Psi(T=\infty)\rangle}$

[TABLE]

III iPEPS representation of thermal states

A projected entangled-pair state state (PEPS) Verstraete and Cirac (2004), also called a tensor product state Nishino et al. (2001); Nishio et al. (2004), is a 2D tensor network representing a state obeying the area law of entanglement. In the simplest case, a PEPS represents a pure state and is built from a network of rank 5 tensors on a square lattice, with one tensor per lattice site. Each tensor has a physical index representing the local Hilbert space of a site. The other four indices of a tensor, called the virtual indices, are contracted with the virtual indices of the neighboring tensors. Their dimension is called the bond dimension $D$ which controls the accuracy of the ansatz. With growing $D$ states with stronger entanglement can be represented by the PEPS. An infinite projected entangled pair state (iPEPS) is a PEPS representing a state on an infinite lattice. To obtain an iPEPS we introduce a unit cell of tensors which is periodically repeated on the lattice, i.e the iPEPS is translational invariant by shifts of the unit cell size. With each site in the unit cell we associate a different PEPS tensor. In this work we use a unit cell with two tensors $A$ and $B$ arranged in a checkerboard pattern (all states studied in this work are compatible with this unit cell).

As $|\Psi(T)\rangle$ is a pure state in the extended Hilbert space, it can be represented by an iPEPS with the tensors having an additional index for the ancillary degrees of freedom, see Fig. 1(a,b). Then $\rho(T)$ (2) can be obtained by a contraction of the iPEPS representation of $|\Psi(T)\rangle$ and its hermitian conjugate, see Fig. 1(c).

The state at infinite temperature, $|\Psi(T=\infty)\rangle$, can be represented exactly by an iPEPS with the bond dimension $D=1$, i.e. a product state (3). For finite $T$ $|\Psi(T)\rangle$ is obtained by an action of an operator $e^{-H/(2T)}$ on the physical degrees of freedom of $|\Psi(T=\infty)\rangle$ (4). As $e^{-H/(2T)}$ in general does not have a numerically tractable exact tensor network representation, we use its Suzuki-Trotter decomposition Trotter (1959); Suzuki (1966, 1976) and VTNR Czarnik and Dziarmaga (2015); Czarnik et al. (2016) to find an iPEPS approximating $|\Psi(T)\rangle$ for a given bond dimension $D$.

IV VTNR

We treat an operator $e^{-H/(2T)}$ as an imaginary time evolution operator $e^{-\tau H}$ with the imaginary time $\tau=1/(2T)$. We decompose $H$ into a sum of classical Hamiltonians $H_{cl}^{j}$, i.e Hamiltonians which are sums of commuting terms,

[TABLE]

We use a second order Suzuki-Trotter decomposition Trotter (1959); Suzuki (1966, 1976) to approximate $e^{-\tau H}$

[TABLE]

The accuracy of the decomposition (5) is controlled by the size of the small time step $\delta\tau$. An exact iPEPO representation of an exponential of a classical Hamiltonian with a finite range of interaction can be found analytically and has $D=d^{2}$ at most (see Ref. Wolf et al., 2008 and simple examples in Refs. Czarnik and Dziarmaga, 2015; Czarnik et al., 2016). Using the Suzuki-Trotter decomposition (5) and iPEPO representations of the exponentials of $H^{j}_{cl}$, we approximate $|\Psi(T)\rangle$ (4) by the 3D tensor network shown in Fig. 2(a).

We use VTNR to approximate the 3D tensor network by an iPEPS with a numerically tractable bond dimension $D$, which yields an approximate representation of $|\Psi(T)\rangle$. In VTNR the iPEPS tensors are obtained by acting with tree tensor networks consisting of isometries Tagliacozzo et al. (2009) on the virtual indices of the 3D tensor network, as shown in Fig. 2(b). The isometries are found by a variational update to minimally distort the partition function $Z(T)$ as described in detail in Refs. Czarnik and Dziarmaga, 2015; Czarnik et al., 2016. The accuracy of the final iPEPS is controlled systematically by the bond dimension of the isometries, which here equals the bond dimension of the final iPEPS.

V Corner Transfer Matrix Renormalization Group (CTMRG)

To compute expectation values of observables and the correlation length $\xi$ we use CTMRG Baxter (1978); Nishino and Okunishi (1996); Orús and Vidal (2009); Corboz et al. (2014). CTMRG approximates contractions of an infinite number of copies of the PEPS tensors by contractions of a finite number of environment tensors $C,E$, where the accuracy is systematically controlled by the bond dimension $\chi$ of the environment tensors. An example of such an approximation is shown in Fig. 3(a,b,c). Details of the algorithm can be found in Refs. Corboz et al., 2014; Czarnik et al., 2016.

VI $\xi$ extrapolation

The correlation length $\xi$ converges very slowly with increasing $\chi$, unlike local observables, e.g. the energy or the magnetization Rams et al. (2018); Corboz et al. (2018). Therefore, to determine the correlation length $\xi$ of $|\psi(T)\rangle$ we use an extrapolation procedure from Ref. Rams et al., 2018 which we summarize briefly below in the simplest case of a translationally invariant PEPS. The procedure uses the eigenvalues of the CTMRG transfer matrix $\mathcal{T}$, see Fig. 3(d).

To set up the extrapolation we define $\epsilon_{j}$ for each eigenvalue $\lambda_{j}$ of $\mathcal{T}$,

[TABLE]

where $j$ numbers the eigenvalues of $\mathcal{T}$ ordered by the absolute value $|\lambda_{0}|\geq|\lambda_{1}|\geq|\lambda_{2}|\dots$ and $\phi_{j}\in(-\pi,\pi]$ determines the phase of $\lambda_{j}$. A connected correlation function of an one-site operator $o$ at a distance $R$ is then expressed as

[TABLE]

Here the form factors $f_{j}^{oo}$ are defined as

[TABLE]

$\mathcal{T}_{o}$ is a transfer matrix of the operator $o$ shown in Fig. 3(e), and $|\phi_{j})$, $(\phi_{j}|$ are left and right eigenvectors of $\mathcal{T}$ normalized as $(\phi_{i}|\phi_{j})=\delta_{ij}$. Therefore the correlation length obtained from $\mathcal{T}$ equals

[TABLE]

As the spectrum of $\mathcal{T}$ is continuous in the limit of $\chi\to\infty$ the extrapolation uses its deviation from continuity $\delta$ as a measure of finite $\chi$ effects. To extrapolate we use eigenvalues of $\mathcal{T}$ contributing to a connected correlation function of the phase transition’s order parameter $C^{mm}(R)$, i. e. the eigenvalues with non-zero form factors of the order parameter $f_{j}^{mm}$. We denote $\epsilon$’s of such eigenvalues by $\{\epsilon^{mm}_{k},k=1,2,\dots\}$ with $\epsilon_{1}^{mm}\leq\epsilon_{2}^{mm}\leq\dots$. We remark that in the case of a second order phase transition the diverging $\xi$ is associated with the symmetry breaking, so we expect that the leading eigenvalue of the transfer matrix determines the asymptotics of the order parameter correlation function, i.e. $\epsilon_{1}^{mm}$ = $\epsilon_{1}$. We observe that this is indeed the case for the transitions investigated below. We note that $\lambda_{0}$ does not contribute to $C^{mm}(R)$ by definition. We define $\delta$ as a distance in between two dominant $\epsilon^{mm}_{k}$, i. e.

[TABLE]

We note that this choice of $\delta$ was proposed and benchmarked in Ref. Rams et al., 2018. Using the spectra from different values of $\chi$, we extrapolate $\epsilon^{mm}_{1}$ as a function of $\delta$ by fitting

[TABLE]

where $\epsilon_{1}^{e}$ is an extrapolated value of $\epsilon^{mm}_{1}$ and $a,b$ are parameters of the fit. The extrapolation gives us

[TABLE]

In Fig. 4 we present examples of the extrapolation for the Quantum Ising model with parameters which are investigated later in Sec. IX.

VII Finite correlation length scaling

Finite Correlation Length Scaling (FCLS) was introduced for simulations of 1D quantum critical phenomena Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012) with infinite MPS. In this case a finite iMPS bond dimension $D$ introduces a finite correlation length $\xi_{D}$ at a critical point. $\xi_{D}$ acts as a cutoff on the diverging correlation length, similarly as a finite system size. It was shown that a scaling analysis in $\xi_{D}$ can be done in a similar way as in conventional finite size scaling, by replacing the system size $L$ by $\xi_{D}$ Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012) in a scaling ansatz, and then make use of this ansatz to obtain the location of the critical point and the values of universal critical exponents. As the finite $D$ introduces also a finite entanglement entropy at the critical point Vidal et al. (2003) FCLS for 1D critical phenomena was originally called Finite Entanglement Scaling Tagliacozzo et al. (2008); Pollmann et al. (2009). The FCLS was recently applied to iPEPS simulations of Lorentz-invariant quantum critical points Corboz et al. (2018); Rader and Läuchli (2018) for which it was found that a finite $D$ introduces finite $\xi_{D}$ at the critical point Corboz et al. (2018).

Here we consider a second order finite temperature phase transition for a quantum Hamiltonian. We use VTNR to find finite $D$ iPEPS approximating purifications of thermal states in the vicinity of the critical temperature $T_{c}$. We observe that the VTNR optimization introduces a finite correlation length $\xi_{D}$ at $T=T_{c}$ (or equivalently at $t=(T-T_{c})/T_{c}=0$) for all $D$’s reached in this work

[TABLE]

This observation motivates us to consider an application of FCLS to obtain $T_{c}$ and the critical exponents from the VTNR results.

We obtain the FCLS ansatz for observables from the standard Finite Size Scaling ansatz by replacing the finite system size by $\xi_{D}$, e. g. for the order parameter $m(t,D)$ we use

[TABLE]

where $\beta$, $\nu$ are the critical exponents, and $\mathscr{M}$ is a non-universal function.

To compute observables we contract the iPEPS using CTMRG. The finite $\chi$ introduces an effective length-scale $\xi_{\chi}$  Nishino et al. (1996), i.e finite $\chi$ values of the observables are given by an more complicated scaling ansatz depending on both $\xi_{D}$ and $\xi_{\chi}$, e. g.

[TABLE]

where $\bar{\mathscr{M}}$ is another non-universal function. To avoid working with the more complicated ansatz we work in the limit of $\chi\to\infty$ as proposed in Refs. Corboz et al., 2018; Rader and Läuchli, 2018. We observe that there is no need to extrapolate the order parameter in $\chi$ as it converges quickly. On the other hand, to obtain a good estimate of $\xi(t,D)\equiv\xi(t,D,\chi\to\infty)$ we use the $\xi$ extrapolation described in Sec. VI.

VIII $T_{c}$ estimation

Finite Size Scaling usually makes use of the Binder cumulant to locate the critical point without prior knowledge of its critical exponents, but in the case of iPEPS computation of the Binder cumulant is challenging because the 4th-order moment of the order parameter would need to be computed. Instead we apply the $m^{\prime}/m$ collapse introduced in Corboz et al. (2018), which makes use of the derivative of the order parameter $m^{\prime}(t,D)=dm(t,D)/dt$, to find $T_{c}$. A FCLS scaling ansatz for $m^{\prime}$ is

[TABLE]

where $\mathscr{M^{\prime}}$ is a non-universal function Corboz et al. (2018). Using (13) and (15) we obtain

[TABLE]

[TABLE]

where $\bar{\mathscr{P}}$ is another non-universal function. Eqs. (16,17) give us the $m^{\prime}/m$ collapse

[TABLE]

We estimate $T_{c}$ by plotting $y=t\frac{m^{\prime}(t,D)}{m(t,D)}$ versus $x=t\frac{m^{\prime}(t=0,D)}{m(t=0,D)}$ for different choices of $T_{c}$ and finding the one for which $y(x)$ data points obtained with different $D$ collapse best onto a single curve.

IX quantum Ising model - benchmark results

The quantum Ising model is given by the Hamiltonian:

[TABLE]

where $\sigma_{z}$, $\sigma_{x}$ are Pauli matrices. For $\Gamma=0$ the model reduces to the classical Ising model with $T_{c}=2/\textrm{ln}(1+\sqrt{2})\approx 2.269$ and for $\Gamma_{c}=3.04438(2)$ Blöte and Deng (2002) it has a quantum critical point. For $0\leq\Gamma<\Gamma_{c}$ it exhibits a low temperature ferromagnetic phase with the order parameter $m=\langle\sigma_{z}\rangle$, which is separated from the paramagnetic phase by a line of finite temperature second order phase transitions belonging to a 2D classical Ising universality class.

As $\Gamma$ is approaching $\Gamma_{c}$, quantum fluctuations are becoming stronger and $T_{c}$ gets suppressed w.r.t. the classical case of $\Gamma=0$. Therefore we expect that with increasing $\Gamma$ the accurate simulation of the finite temperature transition using tensor networks is becoming more challenging since a larger $D$ is necessary to correctly capture the stronger quantum fluctuations. To examine this more closely we investigate in the following $\Gamma=2.5$ as well as a point close to the quantum critical point $\Gamma=2.9$, for which Quantum Monte Carlo (QMC) Hesselmann and Wessel (2016) gives $T_{c}=1.2737(6)$ and $T_{c}=0.6085(8)$, respectively, corresponding to a reduction in $T_{c}$ with respect to $\Gamma=0$ by a factor of 1.8 and 3.7, respectively.

IX.1 $\Gamma=2.5$

We first consider a case well away from the quantum critical point, $\Gamma=2.5$, to provide a proof of principle of the applicability of FCLS to finite temperature VTNR simulations. In Fig. 5 we present data for the order parameter as a function of temperature in the vicinity of the critical temperature, for bond dimensions $D=5-7$. As expected, we do not obtain a sharp phase transition but we see that the order parameter is systematically reduced with increasing $D$ similarly to the case of finite size effects.

We first attempt to estimate $T_{c}$ by using the known critical exponents of the 2D classical Ising universality class, i.e. $\beta=1/8$ and $\nu=1$. To do that we plot $m\xi_{D}^{\beta/\nu}$ which, according to FCLS, should not depend on $D$ at $T=T_{c}$:

[TABLE]

where $a$ does not depend on $D$. Indeed the $m(t,D)\xi_{D}^{\beta/\nu}$ curves for different $D$’s cross as predicted by FCLS, see Fig. 6. We identify $T_{c}$ as the temperature for which the variance of $m(t,D)\xi_{D}^{\beta/\nu}$ is smallest, obtaining $T_{c}=1.2737(2)$ in agreement with the QMC estimate $T_{c}=1.2737(6)$ Hesselmann and Wessel (2016). The $T_{c}$ uncertainty is obtained by varying the range of $D$ and the range of $\chi$ used for $\xi_{D}$ estimation by extrapolation.

To provide further evidence of FCLS we determine the critical exponents, using $T_{c}=1.2737(2)$ found earlier. We first estimate $\beta/\nu$ using the data obtained at $T=T_{c}$. From Eq. (21) we obtain

[TABLE]

A linear fit to the data on a log-log scale shown in Fig. 7(a) yields $\beta/\nu=0.123(15)$, in agreement with the exact $\beta/\nu=1/8$. The error bar takes into account the $T_{c}$ uncertainty and the statistical error of the fit.

Next, we estimate $\beta$ and $\nu$ by performing data collapses based on data in the vicinity of $T_{c}$, using the scaling ansaetze:

[TABLE]

[TABLE]

Using the ansatz (23) we obtain $\beta=0.127(2)$, $\nu=1.04(4)$, in agreement with the exact universality class, see Fig. 7(b). Ansatz (24) yields $\beta=0.126(5),\nu=1.04(5)$, again in agreement with the exact exponents, see Fig. 7(c). In both cases the uncertainties are obtained by taking into account the $T_{c}$ uncertainty and by varying the data range.

Finally, we show that we can estimate $T_{c}$ and $\beta$ without prior knowledge of the universality class nor the value of $T_{c}$. First, we estimate $T_{c}$ by performing a data collapse using the $m^{\prime}/m$ ansatz (18), which yields $T_{c}=1.273(1)$, see Fig. 8(a). The $T_{c}$ uncertainty is obtained by varying the data range. We estimate $\beta$ by performing a data collapse based on the ansatz

[TABLE]

using for $T_{c}$ the value obtained in the $m^{\prime}/m$ collapse. We obtain $\beta=0.12(1)$, see Fig. 8(b). The uncertainty is obtained taking into account the $T_{c}$ uncertainty and varying the data range.

We remark that while the quality of the obtained results is good, we see some deviations from perfect scaling which may be caused either by corrections to finite size scaling or limitations of VTNR in getting optimal tensors. First we note that VTNR is not guaranteed to give the best iPEPS approximation of the thermal state’s purification for a given $D$, as it does not search directly for the best iPEPS tensor representing thermal state purification. Instead it optimizes a tree tensor network (TTN) of isometries which, applied to virtual indices of the tensor network representing a Suzuki-Trotter decomposition of the purification, gives the iPEPS approximating the purification. While this approach makes the variational optimization of the iPEPS efficient it is not equivalent to the most general iPEPS variational optimization procedure. Still our results demonstrate that the accuracy of the optimized tensors is high enough to extract the critical coupling and critical exponents with a good accuracy.

Second we note that for $\Gamma=2.5$ CTMRG convergence is challenging, as for $D=5-7$ we obtain $\xi_{D=7}\sim 350$. For an iPEPS with such a large $\xi$ many iterations of the CTMRG procedure are necessary to converge $m$. Good convergence of the CTMRG environment is important for the variational optimization since VTNR uses the CTMRG environmental tensors to find the best iPEPS 111To ensure a good CTMRG convergence we require the change of $m$ per CTMRG iteration to be smaller than $10^{-8}$. Here to optimize the iPEPS we use $\chi=8D$. We check that using $\chi=6D$ and $\chi=7D$ we obtain results ($T_{c}$ and $\beta$ obtained by the $m^{\prime}/m$ collapse and collapse (25)) in agreement with the ones obtained with $\chi=8D$. Nevertheless we cannot fully exclude the possibility that finite $\chi$ effects contribute to the observed small deviations from the perfect collapse as simulations with larger $\chi$ would be computationally very expensive.

To obtain $\xi_{D}$ we contract the final iPEPS with $8D\leq\chi\leq 24D$ and use the extrapolation procedure described in Sec VI. We remark that in the case of VTNR simulations obtaining convergence in the small Trotter time-step $\delta\tau$ is relatively easy as the computational cost of the simulations scales at most as $O(\textrm{log}(1/\delta\tau))$. For $\Gamma=2.5$ we use a second order Trotter decomposition with $\delta\tau=\tau/2^{10}\leq 0.001$, which is small enough to give results converged in $\delta\tau$.

IX.2 $\Gamma=2.9$

For the more challenging $\Gamma=2.9$ case we analyze VTNR results for $D=7-10$, see Fig. 9(a). Using the $m^{\prime}/m$ collapse we obtain $T_{c}=0.609(4)$, see Fig. 9(b). Using this result for $T_{c}$ and performing a collapse with Eq. (25) we obtain $\beta=0.152(8)$, see Fig. 9(c). While the obtained $T_{c}$ estimate agrees with the QMC estimate $T_{c}=0.6085(8)$ Hesselmann and Wessel (2016) the $\beta$ estimate deviates by about $20\%$ from the exact $\beta=1/8$.

Comparing the results obtained with $D=8-10$ and $D=7-10$ suggests that the $\beta$ estimate still depends significantly on the $D$ range as we obtain $\beta=0.145(5)$ for $D=8-10$ 222Contrary to the case of $\beta$, the $T_{c}$ estimate for $D=8-10$, $T_{c}=0.610(4)$, is similar to the $D=7-10$ estimate. Furthermore, we expect that the necessary $\xi_{D}$ to obtain the asymptotic scaling is larger for $\Gamma=2.9$ than for $\Gamma=2.5$, as $\Gamma=2.9$ is closer to the quantum critical point. Despite larger $D$, the $\xi_{D}$ obtained for $\Gamma=2.9$, although quite large ($\xi_{D}\sim 20-50$), is smaller than $\xi_{D}$ for $\Gamma=2.5$ ($\xi_{D}\sim 40-350$). Therefore we expect that the quality of the results can still be improved by increasing $D$, although it would be computationally very expensive.

The $T_{c}$ and $\beta$ uncertainties are estimated in the same way as in the $\Gamma=2.5$ case. Similar values for $T_{c}$ and $\beta$ are obtained with VTNR using an optimization with $\chi=5D,6D,7D$. A second order Trotter decomposition with a time step $\delta\tau=\tau/2^{10}\leq 0.002$ is used.

X Interacting honeycomb fermions - Benchmark results

We consider a model of interacting spinless fermions on a honeycomb lattice Capponi (2017), given by the Hamiltonian,

[TABLE]

Here $c_{i}(c_{i}^{{\dagger}})$ is a fermionic annihilation (creation) operator at site $i$ and $n_{i}=c_{i}^{{\dagger}}c_{i}$ is a fermion number operator. We set $t=1$ in the following. Furthermore, for the purpose of the benchmark we restrict ourselves to the case of half-filling, $n=\sum_{i=1}^{N}\langle n_{i}\rangle/N=1/2$, for which sign-problem free QMC results are available Wang et al. (2014, 2016); Hesselmann and Wessel (2016). The model has a quantum critical point at $V_{c}=1.356(1)$ Wang et al. (2014). For $V>V_{c}$ there is a low temperature phase with a charge density wave (CDW) order, which is separated from a disordered, high temperature, phase by a line of second order finite temperature phase transitions, which belong to the 2D classical Ising universality class Wang et al. (2016). The CDW order parameter is defined as

[TABLE]

where $\langle n_{A}\rangle$ and $\langle n_{B}\rangle$ are the fermion densities on sub-lattices A and B, respectively. In the limit of $V\to\infty$ the model becomes equivalent to the 2D classical antiferromagnetic Ising model. Here we simulate the model for $V=3$ and the more challenging case of $V=2$, which is closer to $V_{c}$.

X.1 $V=3$

For $V=3$ we analyze $D=12,14,16$ VTNR results in the vicinity of $T_{c}$, see Fig. 10(a). Here QMC predicts $T_{c}=0.936(10)$ Hesselmann and Wessel (2016). We determine $T_{c}$ using the $m^{\prime}/m$ collapse obtaining $T_{c}=0.939(4)$ in agreement with QMC, see Fig. 10(b). Furthermore, we obtain $\beta=0.132(8)$ in agreement with the exact $\beta=1/8$ by performing a data collapse using Eq. (25) and by taking $T_{c}$ obtained from the $m^{\prime}/m$ collapse, see Fig. 10(c).

The $T_{c}$ and $\beta$ uncertainties are obtained similarly as for the quantum Ising model. We use $\chi=5D$ to perform the VTNR optimization obtaining results which are consistent with the ones obtained with $\chi=3D$ and $\chi=4D$. We use a second order Suzuki-Trotter decomposition with a time step $\delta\tau=\tau/2^{11}<0.001$.

X.2 $V=2$

Next we analyze $D=12,14,16$ VTNR results for the more challenging $V=2$ case, see Fig. 11(a). Using the $m^{\prime}/m$ collapse we obtain $T_{c}=0.480(5)$ in agreement with the QMC estimate $T_{c}=0.47(1)$ Hesselmann and Wessel (2016), see Fig. 11(b). Using anatz (25) and taking $T_{c}$ found by the $m^{\prime}/m$ collapse, we obtain $\beta=0.108(4)$, see Fig. 11(c). The $\beta$ estimate deviates by about $10\%$ from the exact $\beta=1/8$. We see that the obtained $\xi_{D}$ for $V=2$ ($\xi_{D}\sim 2-4$) is much smaller than $\xi_{D}$ for $V=3$ ($\xi_{D}\sim 10-20$). Furthermore, as $V=2$ is closer to the quantum critical point than $V=3$ we expect that for $V=2$ a larger $\xi_{D}$ is necessary to be in the asymptotic scaling regime and we expect that the accuracy of $\beta$ can be improved by increasing $D$.

The $T_{c}$ and $\beta$ uncertainties are obtained similarly as for the quantum Ising model. Here we use $\chi=5D$ to perform the VTNR optimization. Consistent results are obtained also with $\chi=3D$ and $\chi=4D$. We use a second order Suzuki-Trotter decomposition with a time step $\delta\tau=\tau/2^{8}<0.01$.

XI Conclusions

In this paper we have studied second order finite temperature phase transitions in the 2D quantum Ising (19) and interacting honeycomb fermion (26) models using infinite projected entangled-pair states (iPEPS) to represent thermal states. The iPEPS were obtained by Variational Tensor Network Renormalization (VTNR). We found that at the critical temperature $T_{c}$ the iPEPS correlation length $\xi_{D}$ is finite for the computationally accessible values of the iPEPS bond dimension $D$. Motivated by this observation we investigated the application of Finite Correlation Length Scaling (FCLS) to obtain precise values of $T_{c}$ and universal critical exponents. We found that in the vicinity of $T_{c}$ the order parameter obeys well the expected behavior predicted by FCLS.

The two models studied in this work exhibit second order finite temperature phase transitions for the transverse fields $\Gamma<\Gamma_{c}\approx 3.04438$ and the interaction strengths $V>V_{c}\approx 1.356$, respectively. At $\Gamma_{c}$ and $V_{c}$ second order quantum phase transitions occur at $T=0$. Using FCLS we obtained estimates of $T_{c}$ and the critical exponents in agreement with the QMC results for $\Gamma=2.5$ and $V=3$ which are sufficiently far from the quantum critical points. For $\Gamma$ and $V$ approaching the quantum critical points we observed that the magnitude of $\xi_{D}$ and the accuracy of the critical data become lower for the same values of $D$. Nevertheless we were still able to obtain $T_{c}$ in agreement with the QMC results for the challenging $\Gamma=2.9$ and $V=2$ cases. For these couplings the values of the critical exponents exhibit a dependence on the range of $D$ values used in the scaling analysis, suggesting that larger $D$’s are needed in order to obtain more accurate estimates of the critical exponents.

In summary, our results further demonstrate the usefulness of tensor network simulations for quantum many-body systems at finite temperature, even for the challenging case of a finite temperature continuous phase transition, for which convergence in $D$ can typically not be reached, but which can be systematically studied using FCLS.

Acknowledgements.

We thank Stephan Hesselmann and Stefan Wessel for providing us numerical values of data published in Ref. Hesselmann and Wessel, 2016 and Marek Rams for useful remarks about the manuscript. This research was funded by the National Science Centre (NCN), Poland under project 2016/23/B/ST3/00830 and the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (grant agreement No. 677061).