Cumulants and Scaling Functions of Infinite Matrix Product States

TL;DR

This paper introduces a numerical scheme using cumulants, Binder cumulant, and scaling functions to accurately determine critical points and exponents in quantum phase transitions of infinite matrix product states, demonstrated on various models.

Contribution

The paper presents a novel, efficient method combining cumulants and scaling functions to extract critical exponents from infinite matrix product states in quantum systems.

Findings

Successfully applied to four 1D models showing accurate critical exponents.

Derived a cumulant exponent relation for consistency checks.

Extended the approach to a 2D model on an infinite cylinder.

Abstract

The order parameter cumulants of infinite matrix product ground states are evaluated across a quantum phase transition. A scheme using the Binder cumulant, finite-entanglement scaling and scaling functions to obtain the critical point and exponents of the correlation length and cumulants is presented. Analogous to the scaling relations that relate the exponents of various thermodynamic quantities, a cumulant exponent relation is derived and used to check the consistency and relationship between the cumulant exponents. This scheme gives a numerically economical way of accurately obtaining the critical exponents. Examples of this scheme are shown for four one-dimensional models - the transverse field Ising model, the topological Kondo insulator, the S = 1 Heisenberg chain with single-ion anisotropy and the Bose-Hubbard model. A two-dimensional model is also exemplified in the square lattice transverse field Ising model on an infinite cylinder. These exemplary systems portray a variety of local and string order parameters as well as phase transition classes that can be studied with the scaling functions and infinite matrix product states.

Equations

$$|\psi ⟩=\sum _{\{s_i\}}\left(\dots A^{s_n}{A}^{s_{n+1}}\dots \right)|\dots s_n{s}_{n+1}\dots ⟩$$

$${\kappa }_1=⟨M⟩,$$

$${\kappa }_2=⟨M^2⟩-{⟨M⟩}^2,$$

$${\kappa }_3=⟨M^3⟩-3⟨M^2⟩⟨M⟩+2{⟨M⟩}^3,$$

$${\kappa }_4=⟨M^4⟩-4⟨M^3⟩⟨M⟩-3{⟨M^2⟩}^2$$

$$+12⟨M^2⟩{⟨M⟩}^2-6{⟨M⟩}^4.$$

$$U_4=1-\frac{⟨M^4⟩}{3{⟨M^2⟩}^2}.$$

$$E_i(L+1)=T_{W_{ii}}\left(E_i(L)\right)+\sum _{j<i}{T}_{W_{ji}}\left(E_j(L)\right),$$

$$T_X\left(E(L)\right)=\sum _{s{s}^{\prime }}⟨s^{\prime }|X|s⟩A^{s^{\prime }†}E(L)A^s.$$

$$E_i(L)=\sum _{m=0}^p{E}_{i,m}{L}^m,$$

$$\displaystyle W=\left(\begin{array}[]{cc}I&Z\\ &I\end{array}\right),$$

$$E_1(L+1)$$

$$E_1(L)$$

$$E_{1,0}=T_I(E_{1,0}),$$

$$E_1(L)=E_{1,0}=\tilde{I}.$$

$$E_2(L+1)$$

$$E_2(1)$$

$$E_2(2)$$

$$E_2(3)$$

$$⋮$$

$$E_2(L)$$

$$E_2(L)=E_{2,0}+E_{2,1}L.$$

$$E_{2,0}+E_{2,1}(L+1)$$

$$E_{2,0}+E_{2,1}+E_{2,1}L$$

$$L^0$$

$$L^1$$

$$\sum _{i=1}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩+e_{2,1,1}\tilde{I}$$

$$\sum _{i=1}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩=e_{2,0,1}|{\lambda }_1⟩+\sum _{i=2}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩,$$

$$e_{2,0,1}\tilde{I}+\sum _{i=2}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩$$

$$\sum _{i=2}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩+e_{2,1,1}\tilde{I}$$

$$\text{Tr}\left(\sum _{i=2}^{m^2}{e}_{2,0,i}|{\lambda }_i⟩{\rho }_R\right)$$

$$\text{Tr}\left(e_{2,1,1}\tilde{I}{\rho }_R\right)$$

$$e_{2,1,1}$$

$$⟨M⟩$$

Full text

I Introduction

One of the cornerstones of quantum mechanics is the postulate that all information of a quantum system is contained in its wavefunction. From a practical perspective, storing and manipulating quantum wavefunctions can be a costly ordeal. This is due to the inherent nature of the Hilbert space of quantum systems, specifically, the size of the Hilbert space and its scaling with the number $N$ of particles of $d$ degrees of freedom is $d^{N}$, i.e. it scales exponentially with system size. Thus, capturing the information of merely 300 two-state particles would require $2^{300}$ bits, a number so large that it would bankrupt the visible universe of its information-storage capacity. This motivates the necessity to concisely represent a wavefunction while revealing the important physics in question.

In the field of low-dimensional many-body quantum physics, a class of ansatz known as tensor network states has been successful in faithfully representing the states of various gapped and gapless phases described by Hamiltonians with local interactions. This success is attributed to the entanglement structure that these physical systems possess, namely where the amount of entanglement entropy is constrained by the system’s physical dimension $D$ Eisert . For example, in gapped systems of size $L$, the entanglement entropy $S$ of the ground state is proportional to $L^{D-1}$. This is known as the “area law of entanglement entropy” and it holds true as long as the system possesses a nonzero energy gap. At a critical point however, the system’s energy gap vanishes and the entanglement entropy violates the area law. In 1D, Ref. Vidal has shown that $S$ diverges logarithmically with $L$. This drastic change of behavior in the entanglement entropy and its ease of computation via tensor network states has made it a viable tool in identifying critical points - by continuously tuning a Hamiltonian parameter, the entanglement entropy would gradually increase approaching the critical point and peaks at the critical point.

Unlike a finite system of size $L$, an infinite translationally-invariant system has no notion of boundaries that allows one to specify the size of the system. Thus, the common notion of length in such systems is the correlation length $\xi$ of some quasiparticle excitation. Being inversely proportional to the energy gap, $\xi$ diverges at a critical point. The tensor network state used to represent an infinite translationally-invariant 1D systems is known as an infinite matrix product state (iMPS) given as

[TABLE]

for a one-site unit cell. $A^{s_{i}}$ is an $m\times m$ matrix and the superscript $s_{i}$ represents an element of the $d$-dimensional local Hilbert space at site $i$. The quantity $m$ is called the bond dimension, and it is a controllable parameter that governs the maximum amount of entanglement that the state can possess. Since $m$ is finite, the peak of $S$ at the critical point does not diverge to infinity, but instead saturates to a value that depends on $m$. The parameter value that corresponds to the maximum $S$ is known as a “pseudocritical point” and has been thoroughly studied in Ref. Tagliacozzo where it was shown that $S\propto\log m^{\kappa}$. $\kappa$ is the finite-entanglement scaling exponent - a quantity related to the central charge of the critical point Pollmann . More importantly, by using finite-entanglement scaling, the authors of Ref. Tagliacozzo presented a scheme to extract $\kappa$ from various universal quantities such as the correlation length, entropy, magnetization order parameter, etc. Analogous to finite-size scaling in finite systems, where critical exponents are extracted by scaling data of universal quantities at the critical point with respect to $L$, the scheme describes extracting the critical exponents by scaling data with respect to $m$. In addition to this, the scaling function of the magnetization order parameter relates the magnetization at different values of $m$ to each other across the critical point. This is again analogous to scaling functions in finite systems where such functions relate universal data of different system sizes to each other. Using the known value of the critical point, the authors used this scaling function to determine $\kappa$ by selecting the value of $\kappa$ that gave the best data collapse of the magnetization of several different values of $m$. Though not done in Ref. Tagliacozzo , the scaling functions could also be used to locate and fine-tune both the critical point and $\kappa$ simultaneously. This is motivated by the fact that the magnetization scales as a power law of $m$ at the critical point and thus its ability to detect the critical point is more pronounced than $S$ which scales logarithmically with $m$ at the critical point. More importantly, generalizing and applying these scaling functions to other universal quantities allows one to determine other critical exponents such as the correlation length exponent $\nu$ and higher order cumulant exponents - all of which are presented in this work.

While an order parameter $M$ is an invaluable tool in discerning phases and locating critical points, there is more information available in its distribution function. Even though this distribution is typically difficult to obtain, it is still possible to gain information about the distribution from the cumulants $\kappa_{n}$ and the higher moments $\mu_{n}=\braket{M^{n}}$ of the order parameter. In probability theory, the cumulants specify the shape of a given distribution. The first cumulant is the mean value

[TABLE]

the second cumulant is the variance

[TABLE]

the third cumulant is the skewness

[TABLE]

and the fourth cumulant is the kurtosis

[TABLE]

Besides probability theory, a modification of the fourth cumulant has found practical usage in the field of phase transitions in what has been known as the Binder cumulant Binder

[TABLE]

By tabulating data the Binder cumulant for different system sizes $L$ across a phase transition, the critical point is read off the point where the Binder cumulant of different system sizes cross each other. The benefit of using the Binder cumulant is that by finding the crossing point between different systems sizes and using two higher order cumulants simultaneously, finite-size effects are drastically reduced. Hence, the critical point obtained from it is much more precise than using solely the order parameter. The occurrence of a crossing point in the Binder cumulant is also observed in infinite translationally-invariant systems represented by an iMPS, where the crossing occurs between $U_{4}$ of different bond dimension $m$ instead of $L$ West .

A true phase transition only occurs when a system is in the thermodynamic limit $L\rightarrow\infty$ where the correlation length $\xi$ diverges. In a finite-sized system, $\xi$ is upper-bounded by the system size $L$. The conventional method of studying a phase transition with a finite MPS is to obtain data in the vicinity of the critical point for several system sizes. Since $\xi$ also depends on $m$, $m$ must be chosen such that the criteria $\xi\approx L$ is fulfilled for each system size. By employing finite-size scaling, the data is extrapolated with respect to $L$ to obtain the quantity of interest in the thermodynamic limit Nishino . This dependence of $m$ on $L$ complicates the procedure of obtaining data since $m$ has to be obtained to satisfy the condition $\xi\approx L$ for each system size. A simpler and more direct way to probe phase transitions would thus be to use an iMPS. There are two advantages of this over a finite MPS. First, since there is no notion of a system size in an iMPS, one does not have to worry about determining the value of $m$ that satisfies the criteria $\xi\approx L$. Hence one only has to extrapolate data with respect to $m$ in order to obtain the data in the $m\rightarrow\infty$ limit. Second, the absence of boundaries completely removes any Friedel oscillations that affects quantities that are dependent on spatial properties of the system.

In the spirit of determining the critical point and exponents via finite-entanglement scaling and scaling functions of the order parameter, this work extends the scheme of Ref. Tagliacozzo through the addition of the order parameter cumulants, the Binder cumulant, and the cumulant scaling functions. The advantage this has over the previous schemes is that this scheme requires a much smaller $m$ in order to determine the critical point and exponents with a higher accuracy. To this end, several 1D examples are given, namely the transverse field Ising model, the topological Kondo insulator, the $S=1$ Heisenberg chain with single-ion anisotropy and the Bose-Hubbard model. A 2D square lattice transverse field Ising model on an infinite cylinder is also investigated. These examples show the capacity of the scaling functions of higher order cumulants in determining the critical point and exponents of a variety of order parameters and phase transitions classes.

II Higher-order Moments and Cumulants in Tensor Networks

To obtain moments and cumulants of an observable, one has to compute the expectation value of the observable operator to different orders. This can be done efficiently through the use of a triangular matrix-product operator (MPO) and its fixed-point equation. The latter is a recursive formula that relates the different elements of the environment matrix $E$ for a given triangular MPO. The advantage this recursive method has over other tensor network methods of calculating higher-order moments and cumulants such as that in Ref. West is that it unifies the method of computing the expectation values of local and string operators whereby the latter ends up being treated as a second order operator (i.e. a two-point correlator). For an upper triangular MPO, this recursive formula is given by Michel

[TABLE]

where $E_{i}$ is the $i$th element of $E$, $L$ are the number of sites, and $W_{ii}$ ($W_{ji}$) are the diagonal (off-diagonal) elements of the triangular MPO. For a lower triangular MPO, the index in the second sum is simply swapped, i.e. $j>i$. $T_{X}$ is the transfer operator that acts the matrix-valued $X$ operator element of the MPO:

[TABLE]

Eq. 7 specifies the action of adding one site to the expectation value in terms of the polynomial form for the $d$ different matrices $E_{i}$ ($1\leq i\leq d$), for a $d\times d$ dimensional MPO McCulloch1 ; McCulloch2 ; Michel . A triangular MPO with zero momentum is characterized by diagonal elements $W_{ii}$ that are proportional to the identity operator, $W_{ii}=xI$, with prefactor $x$ satisfying either $x=1$ or $|x|<1$. As such, the expectation value is a polynomial function of $L$ with matrix-valued coefficients Michel :

[TABLE]

where $p$ is the polynomial degree of $E_{i}(L)$.

The expectation value of operator $M$ written in the form of an upper triangular MPO is $\braket{M}=\text{Tr}\left(E_{d}(L)\rho_{R}\right)$, where $\rho_{R}$ is the reduced density-matrix of the right bipartition of the state, and $E_{d}(L)$ is the component of the left environment matrix that contains the matrix-valued operator that represents the accumulation of all terms of the MPO, corresponding to the $d^{\text{th}}$ column of an upper triangular MPO. For a lower triangular MPO, the expectation value is $\braket{M}={Tr}\left(E_{1}(L)\rho_{R}\right)$ where $E_{1}(L)$ is the component of the left environment matrix that contains the matrix-valued operator accumulated from the first column of the lower triangular MPO. The choice of left (right) environment matrix, left (right) transfer operator and right (left) reduced density matrix is arbitrary, and one could use either choice. As an example, the magnetization order parameter is shown here. This order parameter is given by the upper triangular MPO

[TABLE]

where $Z\equiv S^{z}$ is the $z$ component of the spin operator, and $I$ is the identity. For this triangular MPO, the operator representing the observable is $E_{2}(L)$ which has the expectation value of $\braket{M}=\text{Tr}(E_{2}(L)\rho_{R})$. This quantity is actually trivial to compute via conventional tensor network methods since it is a local expectation value and the MPO is a single-site operator. However, its evaluation from first principles is shown here since, in contrast to the conventional tensor network approach, this method is applicable to more complicated MPO’s such as higher-order moments and string order parameters. The goal here is to show that the correct polynomial degree matrix-valued coefficient Eq. 9 is related to $\text{Tr}(E_{2}(L)\rho_{R})$ – this is done as follows. $W$ has dimension $d=2$, so Eq. 7 gives two terms. First,

[TABLE]

The only operator acting here is transfer operator $T_{I}$ containing the identity operator $I$. As a result, its operation on $E_{1}(L)$ is trivial and therefore independent of $L$. Hence, the polynomial degree $p$ in Eq. 9 for $E_{1}(L)$ is $p=0$ and thus $E_{1}(L)$ can be written as

[TABLE]

Inserting this into Eq. 13 gives

[TABLE]

which implies that $E_{1,0}$ is an eigenvector of operator $T_{I}$ with eigenvalue equal to 1. If the iMPS is in left-canonical form, then $E_{1,0}\propto\tilde{I}$ where $\tilde{I}$ is the $m\times m$ identity matrix. An obvious choice of proportionality factor is to fix $\text{Tr}(E_{1,0}\rho_{R})=1$, which implies

[TABLE]

The second term of the recursion formula is

[TABLE]

Using Eq. 16 and the elements of the MPO from Eq. 12, this becomes

[TABLE]

where in the last step, $C_{Z}\equiv T_{Z}(\tilde{I})$ is a constant matrix, i.e. it has no dependence on any $E_{i}$’s and thus its value can be calculated beforehand. With this, one can now proceed to find an ansatz for the form of $E_{2}(L)$ in the asymptotic large-$L$ limit. This limit shouldn’t depend on the “boundary” form $E_{2}(0)$, hence one can choose $E_{2}(0)=0$ (this choice has no effect on the final solution, up to an irrelevant constant), which leads to

[TABLE]

where $T^{n}_{I}(C_{Z})$ means $n$ repeated applications of the transfer operator, $T_{I}(T_{I}(\ldots(C_{Z})\ldots))$. Thus $E_{2}(L)$ is the form of a geometric series, and the large $L$ limit depends on the nature of the spectral decomposition of $C_{Z}$. Since $T_{I}$ only has a single eigenvalue equal to 1 and all other eigenvalues are strictly less than 1, it follows that $E_{2}(L)$ can diverge at most linearly with $L$, hence the polynomial degree $p$ of Eq. 9 for $E_{2}(L)$ is $p=1$ and thus $E_{2}(L)$ can be written as

[TABLE]

Inserting this into Eq. 18 gives

[TABLE]

Equating powers of $L$ gives

[TABLE]

Similar to Eq. 15, Eq. 23 implies $E_{2,1}\propto\tilde{I}=e_{2,1,1}\tilde{I}$, where $e_{2,1,1}$ is a proportionality constant. The first two indices in the subscript of $e_{2,1,1}$ denotes that it corresponds to $E_{2,1}$ (i.e. same subscript indices for clarity), while the third subscript denotes the eigenvalue number of the transfer operator. Since the first eigenvector of the transfer operator $T_{I}$ is $\tilde{I}$, the third index in the subscript of $e_{2,1,1}$ is $1$. The transfer operator $T_{I}$ can be further decomposed into its components through an eigenvalue decomposition $T_{I}=\sum_{i=1}^{m^{2}}\lambda_{i}\ket{\lambda_{i}}\bra{\lambda_{i}}$, where $\ket{\lambda_{i}}$ are the eigenvectors of $T_{I}$. $E_{2,0}$ can also be expanded in terms of the basis $\ket{\lambda_{i}}$, i.e. $E_{2,0}=\sum_{i=1}^{m^{2}}e_{2,0,i}\ket{\lambda_{i}}$. The constant matrix $C_{Z}$ on the other hand is expanded as $C_{Z}\equiv T_{Z}(\tilde{I})=\sum_{i=1}^{m^{2}}C_{Z,i}\ket{\lambda_{i}}$ where $C_{Z,i}$ are elements of the constant matrix $C_{Z}$. Using these in Eq. 22 gives

[TABLE]

where the orthogonality relation $\sum_{i^{\prime}}\left\langle\lambda_{i}|\lambda_{i^{\prime}}\right\rangle=\delta_{ii^{\prime}}$ was used in the last step. By construction, the spectrum of $T_{I}$ contains the eigenvalue $1$ with corresponding eigenvector $\tilde{I}$. Thus, the eigenvalue decomposition of $T_{I}$ and $E_{2,0}$ can be separated into a part parallel to the identity $\tilde{I}$ and remaining parts that are perpendicular to $\tilde{I}$, i.e.

[TABLE]

with $\ket{\lambda_{1}}=\tilde{I}$, and similarly for the decomposition of $\sum_{i=1}^{m^{2}}\lambda_{i}e_{2,0,i}\ket{\lambda_{i}}=e_{2,0,1}\ket{\lambda_{1}}+\sum_{i=2}^{m^{2}}\lambda_{i}e_{2,0,i}\ket{\lambda_{i}}$, with $\lambda_{1}=1$. Eq. LABEL:eqn:L1_2 then becomes

[TABLE]

As stated above, the goal is to show that the correct polynomial degree of Eq. 9 is related to $\text{Tr}(E_{2}(L)\rho_{R})$. In Eq. 20 the polynomial degree $p=1$ is that corresponding to the coefficient $E_{2,1}$. The latter’s proportionality constant is $e_{2,1,1}$ and this can be related to $\text{Tr}(E_{2}(L)\rho_{R})$ via Eq. 26 by multiplying $\rho_{R}$ on the left and right hand sides and taking the trace:

[TABLE]

The spectrum of the left transfer operator $T_{I}$ contains the eigenvalue $\lambda_{1}=1$ corresponding to its left eigenvector $\tilde{I}$ and right eigenvector $\rho_{R}$ (the right transfer operator on the other hand has eigenvalue 1 corresponding to its left eigenvector $\rho_{L}$, which is the reduced density matrix of the left bipartition, and its right eigenvector $\tilde{I}$). As such, the product between the term parallel to the identity and $\rho_{R}$ in Eq. 27 is $\tilde{I}\rho_{R}=1$, while the products of the terms perpendicular to the identity with $\rho_{R}$ give $\ket{\lambda_{i}}\rho_{R}=0$. Eq. 27 thus reduces to

[TABLE]

where from the first to second line of Eq. 28, $e_{2,1,1}$ and $C_{Z,1}$ are just numbers, so their trace are the numbers themselves. Eq. 28 states that the polynomial degree $p=1$, i.e. the degree of $L^{1}$ with coefficient $E_{2,1}$, is related to the desired expectation value $\text{Tr}(E_{2}(L)\rho_{R})$. Though the terms perpendicular to the identity vanished upon multiplication of $\rho_{R}$, they can be evaluated separately when they are needed for the evaluation of other terms such as in the case of where higher order moments are needed. This derivation can be generalized for the higher order moments and a detailed derivation of $\braket{M^{2}}$ of the MPO in Eq. 12 is shown in Section VIII.1. Ultimately, the moments are expressed as an $n$-degree polynomial in $L$ Michel . For example, the first four moments in terms of cumulants per site are

[TABLE]

Comparing Eq. 29 and the result of Eq. 28 reveals that the coefficient $e_{2,1,1}$ is indeed the first cumulant $\kappa_{1}$.

III Relation between Cumulant Exponents

The fact that the Binder cumulant of different system sizes cross one another at the critical point implies that $U_{4}$ is independent of $L$ at the critical point. This sets a constraint or relationship between the cumulant exponents $\alpha_{i}$ since varying them should not change the value of $U_{4}$ at the critical point. This relationship is revealed by equating the exponent of the bond dimension $m$ in the $n$th order derivative of the singular part of the free energy density $f$ to the exponent of $m$ in the $n$th order cumulant written as a power law function of $m$.

The singular part of the free energy density $f$ introduced in Ref. Privman for a finite system is given by:

[TABLE]

where $L$ is the system size, $d$ is the spatial dimension, $t$ is the reduced temperature given as $t\equiv\frac{T-T_{c}}{T_{c}}$, and $h$ is the scaled applied field $h\equiv H/k_{B}T$. $Y$ is the universal scaling function which even though is universal, depends on system specific properties such as the boundary conditions, lattice geometry, coupling constants, etc, which are captured in the nonuniversal metric factors $C_{1}$ and $C_{2}$. In an infinite system, there is no notion of a system size $L$. Instead, any need for a length is replaced by the correlation length $\xi$, i.e. $L\propto\xi$. This in turn is related to the bond dimension $m$ through the expression $\xi\sim m^{\kappa}$ Tagliacozzo ; Pollmann which will be discussed in the Section IV. One now has $L\propto m^{\kappa}$ and by substituting this into Eq. 33 gives

[TABLE]

From the rules of thermodynamics, the change of the free energy density with respect to a parameter of interest gives the measure of that parameter i.e. its thermodynamic observable or order parameter. Therefore, the first order derivative of $f$ with respect to $h$ gives the first cumulant:

[TABLE]

where from the second to third lines, the Josephson relation $\nu d=2-\alpha$, and the Rushbrooke relation $\alpha+2\beta+\gamma=2$, were combined to eliminate $\alpha$ and give $\nu d=2\beta+\gamma$. The superscript of $Y$ marks the order of the derivative. Proceeding in the similarly fashion for the higher order derivatives of $f$:

[TABLE]

This shows that the exponent of $m$ of any $n$th order cumulant can be obtained from the exponents of just the two first cumulants $\beta$ and $\gamma$.

Cumulants of an order parameter show singular behavior at the critical point, i.e. they either vanish or diverge at the critical point. For an iMPS with a fixed $m$ these divergences still occur, at least in the limit of perfectly converged numerics, although sufficiently close to a pseudo-critical point the iMPS exhibits mean-field-like behavior Liu , i.e. the exponents diverge from their true values. Nevertheless, finite-entanglement scaling of the cumulants at the true critical point follows the correct power-law scaling behavior with the basis size,

[TABLE]

where $\alpha_{i}$ is the cumulant critical exponent. Equating the exponents of $m$ of Eq. 39 and 40 gives a relation between the exponent of the $n$th order cumulant and the first two cumulant exponents:

[TABLE]

Equivalently, a self-contained expression for all $\alpha_{n}$’s in terms of $\alpha_{1}$ and $\alpha_{2}$ can be written as

[TABLE]

where $\alpha_{1}=-\beta\kappa/\nu$ and $\alpha_{2}=\gamma\kappa/\nu$. Eq. 42 is dubbed the cumulant exponent relation and it will be tested for the various models in Section V.

An especially useful application of this relation is in the calculation of certain cumulant exponents whose cumulants are not directly accessible, from those that are readily obtained. For example, in the case when certain symmetries are enforced on a wave function, the odd $n$ order cumulants may turn out to be identically zero, whereas the even $n$ order cumulants are nonzero. In such a circumstance, Eq. 42 can be used to calculate the odd order cumulant exponents from the even ones. To illustrate this, supposing $\alpha_{2}$ and $\alpha_{4}$ have be determined, then by setting $n=4$, Eq. 42 gives

[TABLE]

which gives the first cumulant exponent:

[TABLE]

Subsequently, by setting $n=3$, Eq. 42 gives

[TABLE]

which upon substitution of $\alpha_{1}$ above gives

[TABLE]

IV Finite-entanglement scaling and scaling functions

The idea of finite-entanglement scaling (FES) for infinite systems represented by iMPS was directly adopted from finite-size scaling (FSS) of finite systems Tagliacozzo . The need for FES can be appreciated by understanding the role $m$ plays in the iMPS. To do so, one has to apply a Schmidt decomposition of a state $\ket{\psi}$ which will allow one to observe the entanglement within the bipartition of the state through means of the Schmidt values $\lambda_{n}$. Splitting a state into two parts $A$ and $B$, the Schmidt decomposition is a basis choice that expresses $\ket{\psi}$ as a sum of product between the two parts of the wave function:

[TABLE]

where $\ket{\phi_{n}^{A}}$ and $\ket{\phi_{n}^{B}}$ are orthonormal bases of the Hilbert spaces of the two subsystems $A$ and $B$. The entanglement entropy $S$ between the two subsystems is given as

[TABLE]

Since the size of each matrix $A^{s_{i}}$ in an iMPS is bounded by the bond dimension $m$, all properties of the state at and away from a critical point is dictated by $m$. Away from the critical point, the entanglement entropy $S$ of a ground state is finite, therefore a finite $m$ iMPS would be able to represents the state well. At the critical point, the maximum $S$ for $m$ Schmidt eigenvalues is $\log(m)$. As a result, the finite $m$ dimensional matrices $A^{s_{i}}$ cannot approximate the ground state well and all singular behavior of observables or thermodynamic quantities are blurred out. This is further corroborated by the fact that all correlation functions of MPS decay exponentially, implying that they have finite correlation lengths Ostlund ; Rommer . In spite of this, it is possible to quantitatively measure how observables and thermodynamic quantities behave at and in the vicinity of the critical point through means of FES.

The nature of the transition in a finite and infinite system can be appreciated from their energy landscapes and magnetization order parameter. In a finite system, the energy landscape contains minimums which cross each other at the $m$ and $L$ dependent pseudocritical point $B_{c}(m,L)$. This causes a discontinuous change of the magnetization at $B_{c}(m,L)$ which is a first order transition. As $L$ is increased, these energy minimums move closer to each other in parameter space and finally coincide at $B_{c}(m,L\rightarrow\infty)$ when $L\rightarrow\infty$. This causes the magnetization to vanish smoothly which is a continuous (second-order) phase transition. In an infinite system, the behavior of the order parameter in the vicinity of the critical point is always mean-field in nature. Increasing $m$ shrinks this mean-field behavior around the critical point and the true transition type is recovered in the region left behind by the mean-field behavior. In the case of the transverse field Ising (TFI) model, Ref. Liu has shown this change of the magnetization exponent which is $\beta=1/2$ when $m$ is small whereas a region of $\beta=1/8$ increases in size as $m$ is increased. The former is the well-known magnetization exponent of the TFI model under a mean-field treatment whereas the latter is the true magnetization exponent of the Ising class.

In a finite system, all observables do not form singularities. Taking the correlation length as an example, $\xi$ cannot exceed the system size $L$. Thus, as $\xi$ approaches its maximum value of $\sim L$, it forms a smooth, rounded peak at the $L$-dependent pseudocritical point. This pseudocritical point is located some distance from the true critical point and approaches it as $L$ is increased. This smoothness ensures that $\xi$ is always continuous i.e. no singularity. In contrast to this, the infinite system’s observables always forms a singularity because there is no system-size restriction. Hence, at its $m$-dependent pseudocritical point, $\xi$ forms a sharp, discontinuous peak whose height and pseudocritical point is dependent on $m$. Using these facts, one is able to determine the true critical point by tabulating the singular part of the observable with respect to $m$ and extrapolating the observable or parameter value to the $m\rightarrow\infty$ limit in order to obtain the true value of that observable and the true critical point.

To draw parallels between FSS and FES, let’s first look at FSS. In the thermodynamic limit, the correlation length $\xi(B)$ in the vicinity of the critical point $B_{c}$ diverges as a power law function

[TABLE]

where $\nu$ is the critical exponent of $\xi$. Similarly, suppose there is some observable $\kappa_{n}(B)$ that also diverges at the critical point $B_{c}$ as a power law function:

[TABLE]

where $q$ is the critical exponent of $\kappa_{n}$. In a finite system of size $L$, both $\xi$ and $\kappa_{n}$ are now a functions of $B$ and $L$ i.e. $\xi(L,B)$ and $\kappa_{n}(L,B)$, and they both form peaks at the pseudocritical point $B^{*}(L)$. Since the system is finite, the ground state is effectively ordered when $\xi(L,B^{*}(L))\approx L$ since $\xi(L,B^{*}(L))$ cannot grow larger than $L$. Thus Eqs. 45 for a finite system is written as

[TABLE]

where $g$ is the proportionality constant. Eq. 47 states that the effective distance between the pseudocritical point $B^{*}(L)$ and the true critical point $B_{c}$ is determined by the system size $L$. This equation can be re-expressed in two useful ways. First, as

[TABLE]

where $g^{\prime}=g^{1/\nu}$. This equation states that all $L$-dependent terms on the left hand side must give an overall $L$-independent constant which is on the right hand side. The second expression is

[TABLE]

This is the FSS equation. It enables one to determine $B_{c}$ and $\nu$ from data of an observable e.g. the location of $B^{*}(L)$ corresponding to the peak of $\kappa_{n}(L,B^{*}(L))$, of different values of $L$. Similarly, one can also apply these steps to $\kappa_{n}(B)$ as follows. For a finite system, the peaks of $\kappa_{n}(L,B^{*}(L))$ for a given $L$ occurs $B^{*}(L)$. Using the fact that $\xi(L,B^{*}(L))\approx L$ in the finite system, Eq. 46 becomes

[TABLE]

where $\tilde{g}$ is the proportionality constant. This equation can be re-expressed as

[TABLE]

which states that all the $L$-dependent terms on the left must give an overall $L$-independent constant on the right. This is called the scaling function of $\kappa_{n}(L,B)$. By comparing Eqs. 48 and 51, one notices that the former appears to be a form of scaling of the parameter $B$ while the latter appears to be a scaled form of $\kappa_{n}(L,B)$. Both equations are scaled such that there is no overall dependence of $L$. Therefore by plotting the data of $\kappa_{n}(L,B^{*}(L))/L^{q/\nu}$ versus $|B^{*}(L)-B_{c}|L^{1/\nu}$ for several values of $L$, and by treating the $B_{c}$, $q$ and $\nu$ as tuning parameters, one would find that the data for the different values of $L$ will collapse onto a single curve when the suitable values of the tuning parameters are chosen.

In an infinite system, there is no notion of a system size $L$. Instead, any need of a length is replaced by the correlation length $\xi$. If this system is described by an iMPS, then $\xi$ is related to $m$ by

[TABLE]

where $\kappa$ is the finite-entanglement scaling exponent. This relation was found empirically in Refs. Tagliacozzo ; Andersson , and later derived by Pollmann et al Pollmann . The latter was done by relating the distribution of Schmidt values $\lambda_{n}$ of the bipartition of the wave function, the central charge $c$ of the transition described by the associated conformal field theory, $m$ and $\xi$ at the critical point. This revealed that $\kappa$ is intimately related to the central charge $c$ by:

[TABLE]

The significance of Eq. 53 is that since only a handful of transition classes, and thus $c$’s, are known in 1D, this constraints the possible values of $\kappa$. By using $L\propto\xi\propto m^{\kappa}$, the three important equations 48, 49 and 51 become:

[TABLE]

[TABLE]

and

[TABLE]

Just as in the case of finite $L$, Eqs. 54 and 56 describe a scaled parameter $B$ and scaled function $\kappa_{n}(m,B)$ that are overall independent of $m$ respectively. Hence by plotting data of $\kappa_{n}(m,B^{*}(m))/m^{\kappa q/\nu}$ versus $|B^{*}(m)-B_{c}|m^{\kappa/\nu}$ for several values of $m$, one can expect a data collapse when the suitable values of $B_{c}$, $\gamma$, $\nu$ and $\kappa$ are chosen. Besides the scaling function of the cumulants, Eq. 52 gives a scaling function of $\xi$ that does not have an analog in finite-size scaling. This is written as

[TABLE]

where $g$ is a proportionality constant. Thus plotting this against Eq. 54 would give a data collapse of $\xi(m,B^{*}(m))$ with tuning parameters $B_{c}$, $\nu$ and $\kappa$.

V Results

All ground state wave functions in this work were variationally optimized using the infinite density-matrix renormalization group (iDMRG) algorithm with single-site optimization McCulloch2 ; Schollwock ; Hubig2 . Wave functions of several different bond dimensions $m$ were generated to demonstrate the finite-entanglement scaling of the cumulants $\kappa_{i}$ and correlation length $\xi$, as well as to be used to locate the critical point through means of the Binder cumulant $U_{4}(m)$.

V.1 One-dimensional transverse field Ising model

The 1D transverse field Ising (TFI) model is the quintessential model for studying phase transitions. Its Hamiltonian is given by

[TABLE]

where $B$ is the transverse field strength. The ground state of this model is ferromagnetic when $B<B_{c}$ and paramagnetic when $B>B_{c}$ with $B_{c}=1$. Figure 1 shows the first four cumulants of the magnetization order parameter

[TABLE]

as a function of the transverse field $B$. The ground state is ordered at $B<B_{c}$, giving a nonzero $\braket{M}$, and disordered when $B>B_{c}$, giving $\braket{M}=0$. The variance, skewness and kurtosis on the other hand diverge at the critical point due to large fluctuations in the ground state. One can picture that as the variance diverges, the distribution function expands infinitely. As a result, the skewness and kurtosis of the distribution function also diverges. All four cumulants Figs. 1(a) - 1(d) show the same dependence of $m$, i.e. the point where they vanish/diverge shifts towards the known value of $B_{c}$ as $m$ is increased. This indicates that the iMPS wave function better approximates the true wave function with increasing $m$.

While $B_{c}$ of the 1D TFI model is known, it is good practice to extract its value directly from the cumulant data by employing the Binder cumulant. This also serves as practice for systems where the critical point is not known. In finite systems, the Binder cumulant of several system sizes $U_{4}(L,B)$ are tabulated as a function of a Hamiltonian parameter $B$, and the critical point $B_{c}$ is read off from the point where $U_{4}(L,B)\forall L$ cross or intersect each other, which is denoted $U_{4}(L,B_{c})$. In an infinite system described by an iMPS, the length $L$ is replaced by the correlation length $\xi$ up to a factor $s$:

[TABLE]

and $\xi$ in turn is related to the bond dimension $m$ through Eq. 52. Thus, just as in the case of finite system sizes, the Binder cumulant of different bond dimensions, denoted $U_{4}(m,B)$, can be tabulated and the critical point read off where $U_{4}(m,B)\forall m$ cross or intersect each other. The factor $s$ is treated as a length scaling parameter, and it affects the Binder cumulant by shifting $U_{4}(m,B)$ at different rates that are dependent on the value of $m$, causing $U_{4}(m,B)$ to cross and/or intersect each other at different values of $B$. The optimal value $s^{*}$ in determining the critical point in this work is defined as the value $s$ that gives the crossing or intersection between $U_{4}(m,B)\forall m$, which is denoted $U_{4}(m,B_{c})$. This will be important where there arises a need to distinguish between $U_{4}(m,B_{c})$ and crossing/intersecting points that are formed between $U_{4}(m,B)$ of several different values of $m$ but not all the different values of $m$ - this will be referred to as “spurious crossing points” and they are disregarded as critical point candidates. Since $s^{*}$ is defined to occur at a crossing/intersection of $U_{4}(m,B)\forall m$, it can be obtained by solving $\frac{\partial U_{4}(m,B)}{\partial m}=0$ for $s$ with a linear solver over the range of $B$ - this is how $s^{*}$ is obtained throughout this work. The Binder cumulant for $s=2$ (top) and $s=s^{*}=5.31$ (bottom) are plotted in Fig. 2, where the latter was obtained using a linear solver. In the top figure, $U_{4}(m,B)\forall m$ appears to intersect at $B=1$. However, upon close inspection as shown in the inset of the top figure, this is not a true intersection because $\frac{\partial U_{4}(m,B)}{\partial m}\neq 0$, instead, it is a small but nonzero number. The fact that this occurs at the known critical point is possibly attributed to the simplicity of the 1D Ising model and further examples of other models will elucidate that this is not generic. Within the region of $1<B<1.007$, there are multiple spurious crossing points - each formed from the crossing between pairs of different $m$’s, and can be disregarded as candidates of the critical point. By gradually increasing $s$, the Binder cumulant gradually shift at different rates. The overall effect is that $U_{4}(m,B)$ in the region of $B<1$ moves upwards while $U_{4}(m)$ in the region of $B>1$ moves downwards as can be seen by comparing the top and bottom figures. This causes the intersection point at $B_{c}$ to gradually change into a crossing point at $B_{c}$, all while maintaining the value of $B_{c}=1$. By tuning $s$ to $s^{*}$, $U_{4}(m,B_{c})$ is achieved at $B=1$ which is taken as the critical point $B_{c}$. Increasing $s$ beyond $s^{*}$ does not shift the crossing point any further, however, it must be once again stressed that this observation of having a stable critical point when $s\neq s^{*}$ is special to the 1D Ising model. The latter will be elucidated in the remaining exemplary systems. The key difference of using the Binder cumulant in an iMPS and a finite MPS is the existence of the scaling parameter $s$ in the iMPS. Though it takes an extra step to determine $s=s^{*}$, this step does not require additional simulation time or computational cost. In fact, the overall simulation time and computational cost is greater in the case of finite MPS since data have to be produced for both $L$ and $m$. The utilization of the linear solver requires a guess input value of $s$. Within a guess range of $s$, the linear solver typically converges to a value of $s^{*}$ that gives the single crossing point of $U_{4}(m,B_{c})\forall m$. Whereas far outside this range, the linear solver will either not converge or produce a value $s^{*}$ that does not give the single crossing point of $U_{4}(m,B)$.

Now that $B_{c}$ has been obtained, the finite-entanglement scaling exponent $\kappa$ and the exponents of the cumulants $\alpha_{i}$ (where $i=1,2,3,4$) have to be extracted before the correlation length exponent $\nu$. By taking the logarithm of Eq. 52, $\kappa$ can be extracted as the gradient of the linear equation $\ln\xi=\kappa\ln m$. This is plotted in Fig. 3(a) where the blue circles are the value of the correlation length at the critical point and the red line is the linear fit whose gradient is $\sim 2.083\pm 0.015$. This is rather close to the previously known value of $\kappa$ which is $\sim 2$ Tagliacozzo . In a similar fashion, the cumulants at $B_{c}$ are also plotted as a function of $m$ to extract their critical exponents $\alpha_{i}$ in Fig. 3(b). The respective critical exponents are $\alpha_{1}=-0.257\pm 0.004$, $\alpha_{2}=1.558\pm 0.015$, $\alpha_{3}=3.37\pm 0.03$ and $\alpha_{4}=5.20\pm 0.04$. Eq. 42 allows one to check the consistency and accuracy of the cumulant exponents $\alpha_{i}$’s. Using the obtained values of $\alpha_{1}$ and $\alpha_{2}$, one gets $\alpha_{3}=-\alpha_{1}+2\alpha_{2}=0.257+2(1.558)=3.373$, which differs from value obtained via linear fit in Fig. 3(b) by $\sim 0.15\%$. This difference stems from the uncertainty in the linear fit used to determine the values of $\alpha_{i}$ in Fig. 3(b). Similarly, for the fourth cumulant, $\alpha_{4}=-2\alpha_{1}+3\alpha_{2}=2(0.257)+3(1.558)=5.188$, which differs from the value obtained via linear fit by $\sim 0.25\%$.

The final exponent left to obtain now is $\nu$. By plotting the cumulants according to the scaling function Eqs. 54 and 56 with the obtained values of $B_{c}$, $\kappa$ and the cumulant exponents $\alpha_{i}$, $\nu$ can be tuned to obtain the best data collapse e.g. by minimizing the cumulants’ sum of residual of the different values of $m$ at $B_{c}$. This is shown in Fig. 4 with the result $\nu=1.000\pm 0.005$.

Using Eqs. 54 and 57, the correlation length $\xi$ can also be used to determine $\nu$. Using the obtained values of $B_{c}$ and $\kappa$, $\nu$ is tuned so that the correlation length’s sum of residual of the different values of $m$ at $B_{c}$ is minimized. This is shown in Fig. 5 with $\nu=1.000\pm 0.005$, which is in agreement with that obtained from the data collapse of the cumulants. Eqs. 54, 56 and 57 allows one to additionally fine-tune $B_{c}$ to obtain a better data collapse. Since the value of $B_{c}=1$ gave a sufficiently good data collapse, no further fine-tuning of $B_{c}$ was needed.

In order to make a connection between the more familiar magnetization order parameter exponent $\beta$, one can compare the term $m^{\alpha_{1}}$ from the cumulant scaling function Eq. 56, to that of the conventional scaling function of the magnetization order parameter used in iMPS, $m^{-\beta\kappa/\nu}$ Tagliacozzo . This comparison implies that $\alpha_{1}=-\beta\kappa/\nu$. Using the obtained values of $\alpha_{1}$, $\kappa$ and $\nu$ gives $\beta=-\alpha_{1}\nu/\kappa=0.1234$, which differs from the known value of $\beta=1/8$ by $\sim 1.3\%$. Proceeding in the same way, the second cumulant’s exponent $\alpha_{2}$ is related to the exponent of the variance $\gamma^{*}$ by $\alpha_{2}=\gamma^{*}\kappa/\nu$. Using the obtained values of $\alpha_{2}$, $\kappa$ and $\nu$ gives $\gamma^{*}=\alpha_{2}\nu/\kappa=0.748\sim 3/4$. It is important to note that the second cumulant here is the variance of the order parameter but it has no direct relation to the susceptibility as in the case of the 2D classical Ising model in a longitudinal field. This is because there is no quantum analog of the fluctuation-dissipation theorem that relates the variance to the susceptibility as in the classical case. In the 2D classical Ising model, the susceptibility/variance exponent is $\gamma=7/4$ which can be related to the 1D quantum Ising model with an applied longitudinal field. The exponent $\gamma^{*}=3/4$ obtained in this work is consistent with that found in Ref. Pfeuty .

V.2 One-dimensional topological Kondo insulator

The 1D topological Kondo insulator (TKI) is an effective model introduced in Ref. Alexandrov to study the effects of the strong electron interaction on the topological properties of a 3D bulk insulator. The 1D model consists of a Hubbard chain coupled to a spin-$\frac{1}{2}$ Heisenberg chain by a nonlocal coupling as shown in Fig. 6. The Hamiltonian of this system is given by

[TABLE]

where

[TABLE]

is the 1D Hubbard Hamiltonian (top chain in Fig. 6) describing fermions hopping between sites $j$ and $j+1$ with amplitude $t$, and a Hubbard repulsion $U$ between fermions of opposite spins at site $j$. The second term

[TABLE]

is the 1D Heisenberg Hamiltonian (bottom chain in Fig. 6) that describes the spin exchange between nearest-neighbour $S=\frac{1}{2}$ localized spins. The third term $H_{K}$ represents a non-local Kondo coupling between the Heisenberg and Hubbard chains:

[TABLE]

where $S^{\pm}_{j}$ and $\pi^{\pm}_{j}$ ($S^{z}_{j}$ and $\pi^{z}_{j}$) are the ladder operators ($z$ components) of the spin $\vec{S}_{j}$ in the Heisenberg chain and the p-wave spin density $\vec{\pi}_{j}$ in the Hubbard chain. The $\vec{\pi}_{j}$ operator is given as

[TABLE]

where $\vec{\sigma}$ is the vector of Pauli matrices, and

[TABLE]

The last term in Eq. (61) is the conventional Kondo coupling between a fermion and a localized spin at site $j$,

[TABLE]

In this work, the system is set at half-filling and Hamiltonian parameters that are held fixed are $U=0$, $t=J_{H}=1$, and $J_{K}=2$. The iMPS ground state generated here is enforced with $U(1)$ particle number symmetry to conserve particle number and SU(2) spin rotation symmetry.

When $J_{\perp}<J^{c}_{\perp}$, the ground state is in a symmetry-protected topological (SPT) phase protected by inversion symmetry and undergoes a topological phase transition into a topologically trivial state consisting of local Kondo singlets when $J_{\perp}>J^{c}_{\perp}$ Pillay . This phase transition occurs with the vanishing of the charge excitation gap while the spin excitation gap remains nonzero and finite. This coincides with the presence of a spinless two-particle excitation which is detected with the string order parameter

[TABLE]

where $\hat{n}=\sum_{\sigma}c^{\dagger}_{l,\sigma}c^{l,\sigma}$, which is shown in the Fig. 7(a) as a function of $J_{\perp}$. As shown in Ref. Pillay , $O^{2}_{\text{string}}$ can be expressed as a correlation function $O^{2}_{\text{string}}=\lim_{|j-k|\rightarrow\infty}\braket{p(j)p(k)}$, where $p(j)=\prod_{i<j}(-1)^{\frac{n_{i}-1}{2}}$ is the “kink operator”, and $\braket{p(j)p(k)}=\left\langle\prod_{i<j}\mathbbm{1}_{i}\prod_{i=j}^{k}(-1)^{\frac{n_{i}-1}{2}}\prod_{i>k}\mathbbm{1}_{i}\right\rangle$. $O^{2}_{\text{string}}$ thus appears similar to a local order parameter, e.g., $\braket{M^{2}}$, where $M$ is the magnetization operator. The calculation of $O^{2}_{\text{string}}$ is done by constructing an extensive order parameter $P=\sum_{i}p(i)$, which is written in the matrix product operator (MPO) form as

[TABLE]

The square $\braket{P^{2}}$ is then taken and related to $O^{2}_{\text{string}}$ via $O^{2}_{\text{string}}=\braket{P^{2}}/L^{2}$, where $L=b\xi$, $\xi$ is the correlation length and $b$ is a scaling factor. As shown in Section II, the expectation value of an $n$th power of an operator $P^{n}$ in an iMPS is obtained as a degree $n$ polynomial of the lattice size $L$, which is exact in the asymptotic large-$L$ limit. Thus $O^{2}_{\text{string}}=\braket{P^{2}}/L^{2}$ is evaluated as the coefficient of the degree 2 component of $\braket{P^{2}}$. The example of the evaluation of this second degree polynomial is shown in Section VIII.1.

Fig. 7 shows the cumulants of $O^{2}_{\text{string}}$ as a function of $J_{\perp}$. From these figures, $O^{2}_{\text{string}}$ vanishes, while its variance, skewness and curtosis diverge at $J^{c}_{\perp}\sim 2.21$, consistent with previous results of $J^{c}_{\perp}=2.214$ Pillay . Unlike the case of the 1D TFI model where the magnetization order parameter $\braket{M}$ is finite in the ordered phase and zero in the disordered phase, $O^{2}_{\text{string}}$ remains finite in both SPT and topologically trivial phases and only vanishes at the critical point. This happens because $O^{2}_{\text{string}}$ is not an order parameter in the traditional sense where it discerns an ordered phase from a disordered one. On the contrary, $O^{2}_{\text{string}}$ was constructed specifically to detect a spinless even number-particle excitation. In the TKI model, this excitation is gapped in both phases but vanishes at the critical point. Compared to the cumulants of the 1D TFI model, the cumulants of the TKI in Fig. 7 appear to only shift vertically and not horizontally with increasing $m$. This occurs because the critical point is located within a range of $J_{\perp}$ that is much narrower than the range of $J_{\perp}$ shown in Fig. 7. One can thus expect to see this horizontal shift if one were to use a finer $J_{\perp}$ grid and zoom-in closer around $J_{\perp}\approx 2.20-2.22$.

Even though $O^{2}_{\text{string}}$ does not vanish in one of the phases being separated by the critical point, it can still be used in the Binder cumulant to locate the critical point. The only drawback to this is that the Binder cumulant formed out of $O^{2}_{\text{string}}$ does not form a crossing point between the different values of $m$ for any $s$. Nonetheless, an intersecting point between the different $m$’s still forms and this suffices in locating the critical point. To illustrate this, $U_{4}(m,J_{\perp})$ is tabulated in Fig. 8 for $s=2$ (top) and $s=s^{*}=5.29$ (bottom), where the latter was obtained using a linear solver. While it is noticeable that the top figure appears to have a crossing point at $J_{\perp}\sim 2.16$, this point is far from the critical point that one would expect from observing the cumulants in Fig. 7. Upon close inspection of the vicinity $J_{\perp}\sim 2.16$, one sees that this is in fact not a good crossing point since not all the different values $m$ cross each other simultaneously as shown in the inset of the top figure. The peak however coincides with the critical point $J_{\perp}^{c}=2.21$, but the separation between the different $m$’s does not make it a good indicator of the critical point. By increasing $s$, this peak gradually vanishes as can be seen by comparing the top and bottom figures. When $s=s^{*}$, as in the bottom figure, $U_{4}(m,J_{\perp})\forall m$ intersect at $J_{\perp}=2.21$, which is taken as the critical point $J_{\perp}^{c}$, in agreement with Ref. Pillay . Since $s=s^{*}$, this intersection point satisfies $\frac{\partial U_{4}(m,B)}{\partial m}=0$, which is unlike the top figure of Fig. 2 of the 1D TFI model where the point at $B=1$ was not a true intersection. The point $J_{\perp}=2.22$ appears to be an intersection point, but upon close inspection, it is not since at that point $\frac{\partial U_{4}(m,B)}{\partial m}\neq 0$ as can be seen in the inset of the bottom figure.

Using the critical point $J_{\perp}^{c}=2.21$ obtained from the Binder cumulant, the critical exponents $\kappa$ and $\alpha_{i}$ are obtained by scaling the correlation length of the charge excitation $\xi$ and the cumulants of $O^{2}_{\text{string}}$ with respect to $m$ as shown in Figs. 9(a) and 9(b) respectively. This gives $\kappa\sim 1.027\pm 0.005$, $\alpha_{1}=-0.1056\pm 0.0005$, $\alpha_{2}=0.701\pm 0.003$, $\alpha_{3}=1.552\pm 0.008$ and $\alpha_{4}=2.46\pm 0.01$. The central charge $c$ of this phase transition has been shown to be $c\sim 1$ Pillay which was obtained through the relation of the entanglement entropy $S$ and $\xi$ given by Calabrese

[TABLE]

The gradient of the linear plot of $S$ versus $\ln\xi$ is thus equal to $c/6$ which can be obtained by a linear fit through the data. Alternatively, Eq. 53 can also be used to determine $\kappa$ from the handful of known values of $c$ and vice versa. Using $c=1$, one obtains $\kappa=1.3441$, which significantly differs from the value obtained from the linear fit in Fig. 9(a) by $24\%$. This demonstrates the difficulty and inaccuracy in determining $c$ from $\kappa$. The cumulant exponents obtained from the linear fits in Fig. 9(b) can once again be checked using Eq. 42. Using the obtained values of $\alpha_{1}$ and $\alpha_{2}$, one gets $\alpha_{3}=-\alpha_{1}+2\alpha_{2}=0.106+2(0.701)=1.508$, which differs from value obtained via linear fit in Fig. 9(b) by $\sim 2.9\%$. Similarly, $\alpha_{4}=-2\alpha_{1}+3\alpha_{2}=2(0.106)+3(0.701)=2.315$, which differs from the value obtained via linear fit by $\sim 5.6\%$. As before, this difference of the exponents stem from the uncertainty of the linear fit used. In addition to that, the larger contributor to this difference is the narrow region of $J_{\perp}$ parameter space that the critical point is located in. This affects the chosen value of the critical point $J_{\perp}^{c}$ to which $\xi$ and $\kappa_{i}$ have their respective exponents extracted from. In other words, one would get a smaller difference between the calculated value cumulant exponent and that obtained directly from the linear fit if a finer grid of $J_{\perp}$ was used in detecting the critical point, or if the critical point $J_{\perp}^{c}=2.214$ obtained in Pillay was used directly. The inability to use a finer grid size to zoom-in closer to the location of $J^{c}_{\perp}$ is due to numerical instabilities that plague the cumulant data when the grid size is smaller than $5\times 10^{-4}$ for this particular model.

Using the exponents $\kappa$ and $\alpha_{i}$ obtained from the linear fits in Fig. 9(a) and 9(b), and the respective cumulant and correlation length scaling functions Eqs. 56 and 57, the value of the correlation length’s critical exponent $\nu$ can be obtained. As before, this is done by tuning $\nu$ such that the respective sum of residual of the cumulants and correlation length of different values of $m$ at $J^{c}_{\perp}$ are minimized. This gives the data collapse plotted in Figs. 10 and 11 with the obtained value $\nu=0.710\pm 0.001$. This value of $\nu$ is within the range of previously obtained values in Ref. Pillay of $\nu^{-}=0.666$ and $\nu^{+}=0.742$, where $\nu^{-}$ ($\nu^{+}$) is the exponent obtained from fitting from below (above) the critical point.

V.3 $S=1$ Heisenberg chain with single-ion anisotropy

The $S=1$ Heisenberg chain is a minimal quantum-magnetic model that demonstrates Haldane’s conjecture of the existence of a spectrum gap in integer-spin antiferromagnetic chains Haldane1 ; Haldane2 . The Hamiltonian of this model is

[TABLE]

where $D$ is the single-ion anisotropy term. This model is known to have several phases in the $(D,J_{z})$ parameter space, namely Néel, Haldane, large-$D$, ferromagnetic and two XY phases Chen1 ; Chen2 ; Degli .

In this work, the Hamiltonian parameters held fixed are $J=J_{z}=1$. This reduces the available phases to three. When $D=0$, the ground state is in the Néel phase which possesses a $Z_{2}$ symmetry. When $D<J$, it is in the Haldane phase with an incomplete $Z_{2}\times Z_{2}$ symmetry. This phase is known to be an SPT phase protected by time-reversal ($\mathcal{T}$) or inversion ($\mathcal{I}$) symmetry (dihedral ($\mathcal{D}_{2}$) symmetry is broken by the single-ion anisotropy). When $D>J$, the ground state is in a topologically trivial large-$D$ phase. The transition point between the Haldane and large-$D$ phases has been determined to high precision to be $D_{c}/J=0.96845(8)$ via finite DMRG with system size up to $L=10000$ and $m=1000$ Hu .

The key distinction between the two nonzero $D$ phases is how their ground states transform under the symmetry operation of the above mentioned symmetries. To appreciate this, one has to first look at how the iMPS transforms under symmetry operations. When an iMPS Eq. 1 is in its canonical form, each local tensor $A_{j}$ can be written as a product of $m\times m$ complex matrices $\Gamma_{j}$ and positive, real, diagonal matrices $\Lambda$ Vidal2 which satisfies the canonical condition

[TABLE]

In the canonical form, the transfer matrix

[TABLE]

has a right eigenvector $\delta_{\beta\beta^{\prime}}$ with eigenvalue 1. Similarly, the transfer matrix

[TABLE]

has a left eigenvector $\delta_{\alpha\alpha^{\prime}}$ with eigenvalue 1. If this iMPS is invariant under a local symmetry $g\in G$ which is represented in the local basis as a unitary matrix $u_{g}$, then the $\Gamma_{j}$ matrices must transform under $\left(u_{g}\right)_{jj^{\prime}}$ such that the product in Eq. 1 does not change (up to a phase). This means that the transformed matrices satisfy Garcia

[TABLE]

where $e^{i\theta_{g}}$ is a phase factor and $U_{g}$ is a unitary matrix that commutes with the $\Lambda$ matrices. $U_{g}$ forms an $m$-dimensional projective representation of the symmetry group

[TABLE]

where $\rho\left(g,h\right)$ is the factor set of the representation which can be used to differentiate an SPT phase from a trivial one. Taking time-reversal symmetry as an example, $u_{\mathcal{T}}=U_{\mathcal{T}}K$, where $U_{\mathcal{T}}$ is a basis-dependent unitary e.g. $U_{\mathcal{T}}=e^{i\pi S^{y}}$ for the spin basis, and $K$ is the complex conjugation operation. As such, $\Gamma_{j}$ transforms as

[TABLE]

Relating this to the canonical condition Eq. 74, one finds

[TABLE]

Thus $U_{\mathcal{T}}U_{\mathcal{T}}^{*}$ is an eigenvector of the transfer matrix $T$ Eq. 75 with eigenvalue 1. Since the only unimodular eigenvalue of $T$ is unity and this eigenvalue is unique, one gets $U_{\mathcal{T}}U_{\mathcal{T}}^{*}=e^{i\phi_{\mathcal{T}}}$. Using the unitary property $\mathbbm{1}=U_{\mathcal{T}}U_{\mathcal{T}}^{\dagger}$ and its transpose $\mathbbm{1}=\left(U_{\mathcal{T}}U_{\mathcal{T}}^{\dagger}\right)^{T}=U_{\mathcal{T}}^{*}U_{\mathcal{T}}^{T}$, one can eliminate the $U_{\mathcal{T}}$’s in $U_{\mathcal{T}}U_{\mathcal{T}}^{*}=e^{i\phi_{\mathcal{T}}}$ to get $e^{-2i\phi_{\mathcal{T}}}=1$. The latter sets the restriction $\phi_{\mathcal{T}}=0,\pi$. If $\phi_{\mathcal{T}}=\pi$, then $U_{\mathcal{T}}$ is an antisymmetric matrix which causes the entanglement spectrum to be strictly even-fold degenerate i.e. the ground state is in an SPT phase Pollmann2 . Whereas if $\phi_{\mathcal{T}}=0$, $U_{\mathcal{T}}$ is symmetric and there is no restriction on the degeneracy of the entanglement spectrum i.e. the ground state is topologically trivial. $U_{\mathcal{T}}U_{\mathcal{T}}^{*}=\pm 1$ thus acts as a tool to distinguish the SPT phase from the trivial one. However, in order to evaluate the cumulants, one would need an order parameter that differentiates the two phases based on their phase $e^{i\phi_{\mathcal{T}}}$. This can be achieved through a nonlocal order parameter Pollmann3

[TABLE]

The complex conjugation operator $K_{l}$ acts on the local MPS tensor at site $l$ by complex-conjugating it. Unlike $U_{\mathcal{T}}U_{\mathcal{T}}^{*}$ which obtains the phase of the projective representation from the ancillary states of the iMPS, $O^{2}_{\mathcal{T}}$ obtains it from the physical degrees of freedom. Analogous to the evaluation of $O^{2}_{\text{string}}$ Eq. 68, $O^{2}_{\mathcal{T}}$ is evaluated as the coefficient of the degree 2 component of $\braket{P^{2}}$, where $P=\sum_{i}p_{\tau}(i)$ but with the kink operator $p_{\tau}(j)=\prod_{i<j}u_{\mathcal{T}}(i)$. The MPO form of $P$ is given as

[TABLE]

Following the derivation in the Section VIII.1, this results in $O^{2}_{\mathcal{T}}=\text{Tr}\left(U_{\mathcal{T}}\rho_{R}\right)$, where $\rho_{R}$ is the reduced density matrix of the right bipartition. In the limit where $|j-k|\rightarrow\infty$, $U_{\mathcal{T}}$ is the left eigenvector of the generalized transfer matrix corresponding to $\mathcal{T}$:

[TABLE]

with eigenvalue 1. Since $\rho_{R}$ is chosen from the SVD procedure to be a diagonal matrix, and that $U_{\mathcal{T}}$ is antisymmetric in the SPT phase, the diagonal elements of the product $U_{\mathcal{T}}\rho_{R}$ is zero, ergo $O^{2}_{\mathcal{T}}=\text{Tr}\left(U_{\mathcal{T}}\rho_{R}\right)=0$. In the trivial phase, $U_{\mathcal{T}}$ is symmetric. Hence the diagonal elements of $U_{\mathcal{T}}\rho_{R}$ is nonzero, resulting in $O^{2}_{\mathcal{T}}=\text{Tr}\left(U_{\mathcal{T}}\rho_{R}\right)\neq 0$. Fig. 12 shows the first four cumulants of $O^{2}_{\mathcal{T}}$. The first cumulant $\kappa_{1}$ is the nonlocal order parameter $O^{2}_{\mathcal{T}}$ itself which is zero when $D\leq 0.96$, i.e. in the SPT Haldane phase protected by $\mathcal{T}$. On the other hand, $\kappa_{1}=0$ when $D>0.96$ which is the topologically trivial large $D$ phase. The variation with respect to $m$ of $\kappa_{1}$ is minute, but this is more apparent in the other cumulants $\kappa_{2}$ - $\kappa_{4}$, which all diverge at the critical point. Just as in the case of the cumulants $O^{2}_{\text{string}}$ of TKI in Fig. 7, the cumulants of $O^{2}_{\mathcal{T}}$ do not shift horizontally with increasing $m$, indicating that the critical point is located in narrow region of $D$ compared to the range of $D$ shown in Fig. 12. The sharp transition in the cumulants makes identifying the critical point an easy task which is taken to be $D_{c}=0.96$.

Even though the critical point has already been located through the use of the cumulants, the Binder cumulant can still be used as a consistency check of this critical point’s location. Fig. 13 shows $U_{4}$ of $O^{2}_{\mathcal{T}}$ for several values of $s$. As $s$ is increased, $U_{4}$ in the region $D<D_{c}$ ($D>D_{c}$) shifts downwards (upwards), as can be seen in Figs. 13(a) to 13(d). The optimum $s^{*}$ that gives a crossing point $U_{4}(m,D)\forall m$ occurs at $D=0.96$ when $s^{*}=3.12$ as shown in Fig. 13(c). This value of $s^{*}$ was obtained using a linear solver to solve $\frac{\partial U_{4}(m,D)}{\partial m}=0$ for $s$. The critical point is therefore taken to be $D_{c}=0.96$, in agreement with that obtained from the cumulants.

Since $\kappa_{1}$ and $\kappa_{3}$ vanish when $D\leq 0.96$, it is not possible to use them to obtain their critical exponents at $D=D_{c}$. However, it is still possible to use the data in the vicinity of the critical point to obtain the critical exponents. Here, using the data at $D=0.97$, $\xi$ and $\kappa_{i}$ are plotted against $m$ in Fig. 14. By fitting these data with a linear fit and extracting their gradients, the critical exponents are obtained: $\kappa\sim 1.341\pm 0.015$, $\alpha_{1}=-0.133\pm 0.001$, $\alpha_{2}=0.7987\pm 0.005$, $\alpha_{3}=1.737\pm 0.008$ and $\alpha_{4}=2.53\pm 0.02$. From Eq. 53, the closest value of $\kappa$ to the one obtained here that gives a known value of $c$ is $\kappa=1.3441$, corresponding to $c=1$, in good agreement with Refs. Hu ; Chen1 ; Chen2 ; Degli . The cumulant exponents can also be checked using Eq. 42. Using the obtained values of $\alpha_{1}$ and $\alpha_{2}$, one gets $\alpha_{3}=-\alpha_{1}+2\alpha_{2}=0.133+2(0.7987)=1.7304$, which differs from value obtained via linear fit in Fig. 14(b) by $\sim 0.4\%$. Similarly, $\alpha_{4}=-2\alpha_{1}+3\alpha_{2}=2(0.133)+3(0.7987)=2.662$, which differs from the value obtained via linear fit by $\sim 5\%$.

The obtained values of the exponents $\kappa$ and $\alpha_{i}$ are now used in the scaling functions of the cumulants and the correlation length Eqs. 56 and 57 respectively to determine the value of the exponent $\nu$. By tuning $\nu$ such that the sum of residual of the cumulants and correlation length of different values of $m$ at $D_{c}$ is minimized, the best data collapse is obtained which marks the optimum value of $\nu$. For the sake of consistency, the same value $D=0.97$ used to obtain the exponents $\kappa$ and $\alpha_{i}$’s is used here. The data collapse of the cumulants are displayed in Fig. 15 and that of the correlation length is displayed in Fig. 16 where the best data collapse occurs when $\nu=1.470\pm 0.001$. This differs with the value $\nu=1.472(4)$ obtained in Ref. Hu by $\sim 0.16\%$.

V.4 Two-dimensional square lattice transverse field Ising model on an infinite cylinder

The 2D square lattice transverse field Ising model on an infinite cylinder is given by the Hamiltonian

[TABLE]

where $B$ is the transverse field strength. This is a semi-infinite cylinder i.e. its length is infinite but possesses a finite circumference $w$ that the index $j$ sums over. The circumference chosen for this work is 12 sites. Just like in the 1D case, the order parameter is the magnetization $\braket{M}=\sum_{i,j}\sigma^{z}_{i,j}$, which is now summed over two indices $i$ and $j$ since the system is two dimensional. Another similarity shared between the 1D and 2D TFI models is the behavior of $\braket{M}$. When $B<B_{c}$, the ground state is ordered and thus $\braket{M}\neq 0$. Whereas when $B>B_{c}$, the ground state is disordered, hence $\braket{M}=0$. Figure 17(a) - (d) shows the first four cumulants of $\braket{M}$ as a function of $B$. The behavior of the cumulants closely resemble that of the 1D TFI in Fig. 1 where $\kappa_{1}$ vanishes upon approaching the critical point while the other cumulants diverge. The cumulants shift significantly with increasing $m$, where their vanishing/divergence approach the critical point as $m$ is increased.

Since the cumulants do not easily locate the critical point, the Binder cumulant is used to do so. $U_{4}$ as a function of $B$ is shown in Fig. 18 for three different values of $s$. Just as in the previous examples, varying $s$ shifts $U_{4}(m)$ at different rates for the different $m$ above and below the critical point. When $s<s^{*}$, for example in Fig. 18(a) where $s=2$, there are multiple spurious crossing points formed between pairs of different $m$’s. This is analogous to the Binder cumulant for $s=2$ in the top figure of Fig. 2 for the 1D Ising model. Increasing $s$ shifts the values of $U_{4}(m)$ such that $U_{4}(m)$ in the region below the critical point shift upwards while the values above the critical point shifts downwards - again, analogous to the increasing $s$ in the bottom figure of Fig. 2. As a results, nearby spurious crossing points merge to form $U_{4}(m)^{*}$. This can be seen in Fig. 18(b) where $s=s^{*}=3.27$ and $U_{4}(m)^{*}$ occurs at $B=3.01$ which is taken as the critical point $B_{c}$. As before, $s^{*}$ was obtained by solving $\frac{\partial U_{4}(m,B)}{\partial m}=0$ for $s$ using a linear solver. In the region of $3.015<B<3.018$, the $U_{4}(m=500)$ (black inverted triangles) cross the other values of $m$ at, forming multiple spurious crossing points between pairs of $m$’s. When $s$ is further increased beyond $s^{*}$, $U_{4}(m)$ shift at different rates, thus once again forming multiple spurious crossing points between pairs of different $m$’s as can be seen in Fig. 18(c) where the exaggerated value of $s=20$ is chosen to clearly demonstrate this. As before, these spurious crossing points can be disregarded as candidates of the critical point.

Using $B_{c}=3.01$ obtained from the Binder cumulant, the exponents of the correlation length and cumulants are now extracted by plotting $\xi$ and $\kappa_{i}$ against $m$. Fig. 19(a) shows the log-log plot of correlation length vs. $m$ for $B=3.01$. The gradient of this linear fit gives the critical exponents $\kappa=1.024\pm 0.006$. Fig. 19(b), on the other hand, is the log-log plot of the cumulants with their respective linear fits. By comparing the cumulants in this figure, one notices that the higher the cumulant order, the larger the fluctuation with respect to $m$. This makes the error in the linear fit for the higher cumulants larger, and thus the gradient of the higher cumulants is more prone to error. The gradient of the linear fits give the cumulant exponents $\alpha_{1}=-0.276\pm 0.002$, $\alpha_{2}=0.904\pm 0.006$, $\alpha_{3}=2.11\pm 0.01$, $\alpha_{4}=3.110\pm 0.012$. The quality of the linear fit can be checked from the cumulant exponent relation Eq. 42. Using the obtained values of $\alpha_{1}$ and $\alpha_{2}$, one gets $\alpha_{3}=-\alpha_{1}+2\alpha_{2}=0.276+2(0.904)=2.084$, which differs from value obtained via linear fit in Fig. 19(b) by $\sim 1.2\%$. Similarly, $\alpha_{4}=-2\alpha_{1}+3\alpha_{2}=2(0.276)+3(0.904)=3.264$, which differs from the value obtained via linear fit by $\sim 5\%$. As stated earlier, this large difference between the calculated and fitted values of the higher order cumulants’ exponents stem from the fact that the higher order cumulants tend to fluctuate more. As a result, their fitted exponents are more susceptible to errors.

Using the values of $B_{c}$, $\kappa$ and $\alpha_{i}$ obtained, together with the scaling function of the cumulants Eq. 56 and correlation length Eq. 57, the value of $\nu$ can now be obtained. This is done by tuning $\nu$ such that the respective sum of residual of the cumulants and correlation length of the different values of $m$ at $B_{c}$ are minimized. This gives the value $\nu=0.750\pm 0.008$ and the data collapse of the cumulant shown in Fig. 20 and that of the correlation length shown in Fig. 21. This value of $\nu$ sits in between the value of the 1D TFI obtained in Section V.1 where $\nu=1$, and the full 3D classical (or equivalently 2D quantum) Ising model in Ref. Rader where $\nu=0.629970(4)$. This is expected since the infinite cylinder sits geometrically in between the infinite chain (full 1D) and the infinite plane (full 2D).

V.5 One-dimensional Bose-Hubbard model

The 1D Bose-Hubbard (BH) model is given by the Hamiltonian

[TABLE]

where $b^{\dagger}_{j}$ ($b_{j}$) is the boson creation (annihilation) operator at site $j$, $n_{j}=b^{\dagger}_{j}b_{j}$ is the number operator at site $j$, and $U$ is the Hubbard repulsion term that penalizes double occupancy. For one particle per site, and fixing the energy scale $U=1$, previous studies based on a variety of methods such as finite and infinite DMRG, quantum Monte Carlo, exact diagonalization and Bethe ansatz have shown to undergo a phase transition from the Mott insulator to a superfluid phase at $J_{c}\approx 0.26-0.3$ Fisher ; Kuhner ; Gu ; Krutitsky ; Carrasquilla ; Rams (and references therein). The latest work is that of Ref. Rams where $J_{c}=0.3048(3)$ was obtained by using a maximum bond dimension of $m=4000$ and an extrapolation of the correlation length with respect to the first two eigenvalues of the transfer matrix. This transition belongs to the Berezinskii-Kosterlitz-Thouless (BKT) universality class Fisher ; Krutitsky . Unlike other classes of phase transitions where the correlation length diverges algebraically as a power law as the critical point is approached, the BKT transition is known for its exponentially diverging correlation length Kuhner ; Loison

[TABLE]

which causes all data of the critical point to be strongly plagued by finite-size effects.

In this work, the Hubbard repulsion $U$ is set to unity and the system is at half-filling. In addition, $U(1)$ particle number symmetry is enforced so that the total particle number per site is conserved. Even though the Mermin-Wagner theorem forbids the spontaneous-breaking of a continuous symmetry in 1D, an iMPS formulated without explicitly preserving $U(1)$ will spontaneously break this symmetry. For example, consider the $m=1$ limit where the iMPS reproduces the mean-field solution, i.e. the iMPS groundstate is the exact solution of a mean-field Hamiltonian. In this limit, the iMPS explicitly breaks $U(1)$ symmetry. When $J>J_{c}$, this groundstate is comprised of a superfluid where the number of particles at each site is a superposition of all possible particle numbers, i.e. $\ket{\psi}=\left(a_{0}\ket{0}+a_{1}\ket{1}+a_{2}\ket{2}+\ldots\right)^{\otimes L}$. In the large-$m$ limit the $U(1)$ symmetry is restored, and the superfluid order parameter vanishes as required by the Mermin-Wagner theorem. By explicitly preserving $U(1)$ symmetry, the superfluid order parameter is always zero and thus is not a good choice of an order parameter. In contrast, the Mott string order parameter is zero in the superfluid phase but non-zero in the Mott-insulating phase. The latter is the region $J<J_{c}$ where the groundstate is dominated by the Hubbard repulsion term, therefore, double-occupancy is energetically expensive and the groundstate is comprised of a one-particle-per-site insulator. The Mott string order parameter is given as

[TABLE]

where $\hat{n}_{l}=b^{\dagger}_{l}b_{l}$ measures the number of bosons at site $l$. In the Mott insulator phase, the strong repulsion between bosons causes long-range correlations of single-site occupancy (alternatively, fluctuations in the density are short-range correlated), thus $O^{2}_{\text{Mott}}\neq 0$. On the other hand, in the gapless phase, long-range (power-law) fluctuations in the density causes $O^{2}_{\text{Mott}}=0$. Fig. 22(a)-(d) show the first four cumulants of $O^{2}_{\text{Mott}}$. All the cumulants shift towards the critical point as $m$ is increased. As can be seen in Fig. 22(a), $\kappa_{1}$ decreases but does not vanish completely in the gapless superfluid phase. However, $\kappa_{1}$ in this phase decreases with increasing $m$, indicating that it would vanish completely in the $m\rightarrow\infty$ limit. This corresponds to the power-law fluctuations in density being restricted to a finite correlation length by the finite $m$.

The Binder cumulant in finite system studies has been known to be cumbersome in locating the critical point of a BKT transition because one does not simply obtain a single crossing point between the different system sizes Loison . Instead, multiple crossing points between each system size is observed. As a results, one has to compare each crossing point and extrapolate them to get the final crossing point which marks the critical point. Even so, this extrapolation is not straight-forward and can give a very different critical point if one is not careful to do the correct comparison between many different system sizes. Because of this, a large number of system sizes and comparisons are required to locate the correct critical point. Projecting this fact onto the consideration of an infinite system described in this work, the additional length scaling parameter $s$ would add an extra degree of difficulty since there is no clear way to determine its optimum value $s^{*}$ which is defined as the crossing point of the Binder cumulant for all values of $m$. As such, the deployment of the Binder cumulant to determine the critical point is deemed impractical. Instead, the scaling functions of the cumulants and correlation length are used directly to determine the critical point and critical exponents simultaneously.

Since $\xi$ in Eq. 88 scales as an exponential instead of a power law as in Eq. 45, a new form of the scaled parameter Eq. 48 has to be obtained. Starting from Eq. 88 for a finite system and following the same steps to obtain Eq. 48 from Eq. 45, one gets:

[TABLE]

Converting this to an infinite system described by an iMPS by the substitution $L\rightarrow\xi\propto m^{\kappa}$, gives

[TABLE]

where $g$ is a proportionality constant and $g^{\prime}=g^{1/\nu}$. Using this together with the cumulant scaling function of the form Eq. 56 and the correlation length scaling function Eq. 57, the data collapse of the cumulants and correlation length can be obtained respectively. The data collapse of the four cumulants are plotted in Fig. 23. By tuning each cumulant exponents $\alpha_{i}$ simultaneously with the critical point $J_{c}$ and exponents $\kappa$ and $\nu$, the best data collapse, i.e. the collapse corresponding to the cumulants’ minimized sum of residual of the different values of $m$ over a range of $J$, is obtained with the critical point $J_{c}=0.2850\pm 0.0005$, and exponents $\kappa=1.275\pm 0.001$, $\nu=0.500\pm 0.001$, $\alpha_{1}=-0.375\pm 0.001$, $\alpha_{2}=0.350\pm 0.005$, $\alpha_{3}=1.10\pm 0.02$, $\alpha_{4}=1.80\pm 0.05$. Using Eq. 42 to check the consistency of the cumulant exponents, one finds that $\alpha_{3}=-\alpha_{1}+2\alpha_{2}=0.375+2(0.35)=1.075$, which differs by $\sim 2.3\%$ from the value obtained by directly tuning $\alpha_{3}$ in the data collapse scaling function. Similarly, $\alpha_{4}=-2\alpha_{1}+3\alpha_{2}=2(0.375)+3(0.35)=1.80$, which is exactly the value obtained from the data collapse of the scaling function. The critical point obtained here however differs from the value obtained in Refs. Krutitsky ; Carrasquilla ; Rams of $J_{c}\approx 0.3$ by $5\%$. The exponent $\nu$ agrees exactly with that in Ref. Fisher .

The data collapse of the correlation length scaling function Eq. 57 is now used to check the whether the critical point and exponents $\kappa$ and $\nu$ obtained from the cumulant scaling function were correct. By minimizing the correlation length’s sum of residual of different values of $m$ over a range of $J$, the data collapse obtained is shown in Fig. 24 with the values $J_{c}=0.295\pm 0.001$, $\nu=0.500\pm 0.001$, $\kappa=1.275\pm 0.001$. The critical point here is closer to that obtained in Refs. Krutitsky ; Carrasquilla ; Rams of $J_{c}\approx 0.3$, differing by $\sim 1.7\%$. The superfluid phase, comprising of free bosons, is described by central charge $c=1$ Krutitsky . Using this, Eq. 53 gives $\kappa=1.344$ which differs from the obtained value $\kappa=1.275$ by $\sim 5.1\%$.

VI Summary

The order parameter cumulants were studied across the second order and BKT transition classes for several 1D and 2D exemplary systems. These cumulants were obtained using the recursive formula Eq. 7 which offers an efficient way of computing operators of any order and unifies the procedure of calculating both local and string operators. Using the Binder cumulant, finite-entanglement scaling and scaling functions, the critical point and exponents were determined with a relatively smaller bond dimension compared to previously-known data. The procedure to obtain the critical point and exponents are summarized here:

Obtain the first four cumulants of the order parameter and correlation length $\xi$ as a function of Hamiltonian parameter $B$ across the critical point for several values of bond dimension $m$. 2. 2.

Using the four cumulants, compute the Binder cumulant $U_{4}(m,B)$ as a function of $B$ for all values of $m$. 3. 3.

Tune the length scale parameter $s$ relating the system size $L$ to the correlation length $\xi$ in Eq. 60 to obtain the crossing/intersection point of $U_{4}(m,B)$ for all $m$. This crossing point is taken as the critical point $B_{c}$. Alternatively, use a linear solver to solve $\frac{\partial U_{4}(m,B)}{\partial m}=0\,\,\forall m$ over a range of parameter $B$. This will give $B_{c}$ at the optimum value $s=s^{*}$. 4. 4.

Using $B_{c}$, employ finite-entanglement scaling to extract the finite-entanglement scaling exponent $\kappa$ and the cumulant exponents $\alpha_{i}$. This is done by plotting $\log\xi(B_{c})$ and $\log\kappa_{i}(B_{c})$ against $\log m$ respectively. Fit a linear function through these data. The critical exponent is the gradient of the linear fit. Check the consistency of the cumulant exponents by using the cumulant exponent relation Eq. 42. If the exponents disagree significantly, adjust the linear fit parameters. 5. 5.

Use the cumulant scaling function Eq. 56, the correlation length scaling function Eq. 57 and the scaled parameter Eq. 54 to obtain the data collapse of the cumulants and correlation length respectively. Tune the exponent $\nu$ to obtain the best data collapse. Fine-tune the critical point $B_{c}$ if necessary. 6. 6.

If $B_{c}$ is significantly different from that obtained in step 3, go back to step 4, adjust the linear fit parameters, and repeat.

Where the Binder cumulant is cumbersome in producing the critical point, such as in the case of the BKT transition class, the scaling function can be directly employed to obtain the critical point and exponents simultaneously. The consistency of the exponents can be then checked with the cumulant exponent relation.

VII Acknowledgement

I.P.M. acknowledges support from the ARC Future Fellowships Scheme No. FT140100625.

VIII Appendix

VIII.1 Derivation of $\braket{M^{2}}$

This section shows the derivation of the second order moment $\braket{M^{2}}$ for the operator $M$ given in Eq. 12. This is an extension of the procedure used in the calculation of the order parameter shown in Section II and can be easily generalized to any higher order moments.

The operator $M^{2}$ is given by

[TABLE]

As stated in Section II, the expectation value of an operator in the form of an upper triangular MPO is calculated by the fixed-point equation of the environment matrix $E$:

[TABLE]

For the upper triangular MPO, the operator representing the observable is $E_{4}(L)$ which has the expectation value of $\braket{M^{2}}=\text{Tr}(E_{4}(L)\rho_{R})$. The goal here is to show that the correct polynomial degree matrix-valued coefficient Eq. 9:

[TABLE]

is related to $\text{Tr}(T_{Z^{2}}(E_{4}(L))\rho_{R})$ - this is done as follows. Writing out each term $E_{i}(L+1)$,

[TABLE]

The expectation value of interest $\text{Tr}(E_{4}(L)\rho_{R})$ is in Eq. 107. To obtain this, one has to solve the Eqs. 104, 105 and 106 sequentially in order to obtain terms that will be substituted into the final equation Eq. 107.

Starting from $E_{1}(L+1)$, Eq. 104 only contains the operator $T_{I}$, which acts trivially on $E_{1}(L)$ and is therefore independent of $L$. Hence the polynomial expansion of $E_{1}(L)$ is of polynomial order $p=0$:

[TABLE]

Inserting this into Eq. 104 gives

[TABLE]

which implies that $E_{1,0}$ is an eigenvector of $T_{I}$ with eigenvalue unity, and hence $E_{1,0}\propto\tilde{I}$, where $\tilde{I}$ is an $m\times m$ matrix. This result will be used in subsequent equations where $E_{1}$ is needed.

The next equation is Eq. 105 which contains operators $T_{I}$ and $T_{Z}$. Following the same reasoning to obtain the ansatz of $E_{2}(L)$ via Eq. 19 leading to Eq. 20, one finds that the polynomial order of $E_{2}(L)$ is thus $p=1$ and its polynomial expansion is

[TABLE]

Inserting this and $E_{1,0}\propto\tilde{I}$ into Eq. 105 gives

[TABLE]

where $C_{Z}\equiv T_{Z}(\tilde{I})$ is a constant matrix, i.e. it has no dependence on any $E_{i}$’s and thus its value can be calculated beforehand. Equating powers of $L$ gives

[TABLE]

Eq. 113 implies $E_{2,1}\propto\tilde{I}=e_{2,1,1}\tilde{I}$, where $e_{2,1,1}$ is a proportionality constant. Inserting this into Eq. 112 gives

[TABLE]

By decomposing $E_{2,0}=\sum^{m^{2}}_{i=1}e_{2,0,i}\ket{\lambda_{i}}$, $T_{I}=\sum^{m^{2}}_{i=1}\lambda_{i}\ket{\lambda_{i}}\bra{\lambda_{i}}$ and $C_{Z}\equiv T_{Z}(\tilde{I})=\sum^{m^{2}}_{i=1}C_{Z,i}\ket{\lambda_{i}}$, where $C_{Z,i}$ are elements of the constant matrix $C_{Z}$, Eq. 114 becomes

[TABLE]

where the orthogonality relation $\sum_{i^{\prime}}\left\langle\lambda_{i}|\lambda_{i}^{\prime}\right\rangle=\delta_{ii^{\prime}}$ was used in the last step. Further decomposing Eq. LABEL:appendix_eqn:E2_5 into a term parallel to the identity and terms perpendicular to the identity,

[TABLE]

The parts that are parallel to the identity in Eq. 116 are

[TABLE]

Multiplying both sides by of Eq. 117 by $\rho_{R}$ and taking the trace gives

[TABLE]

where $\tilde{I}\rho_{R}=1$ was used. The terms that are perpendicular to the identity in Eq. 116 will be used in the later part for calculating $E_{1}(L)$ where the value of $e_{2,0,i}$ will be needed. Since there is no explicit way of determining $e_{2,0,i}$, it has to be solved with a numerical solver. This is done be rewriting the parts of Eq. 116 that are perpendicular to the identity, i.e.

[TABLE]

as a set of linear equations

[TABLE]

Eq. 120 can now be solved for $e_{2,0,i}$ with a linear solver such as the generalized minimal residual solver (GMRES).

Eq. 106 is similar to that of Eq. 105. Applying the same procedure gives

[TABLE]

and

[TABLE]

This implies

[TABLE]

One can now look for the ansatz for the final term $E_{4}(L)$ in Eq. 107 in the asymptotic large-$L$ limit. To do so, the similar reasoning to obtain the ansatz of $E_{2}(L)$ via Eq. 19 leading to Eq. 20 is applied here to obtain the ansatz for $E_{4}(L)$. In this case, $E_{4}(L)$ is the form of a geometric series, and the large $L$ limit depends on the nature of the spectral decomposition of $C_{Z}$ and $C_{Z^{2}}$. Since $T_{I}$ only has a single eigenvalue equal to 1 and all other eigenvalues are strictly less than 1, $E_{4}(L)$ can diverge at most quadratically with $L$, hence $E_{4}(L)$ has polynomial order $p=2$ and its polynomial expansion is

[TABLE]

Inserting this, together with $E_{1}(L)$, $E_{2}(L)$ and $E_{3}(L)$, into Eq. 107 gives

[TABLE]

Equating powers of $L$,

[TABLE]

where the last term on the right hand side of Eq. 126 was $T_{Z^{2}}(E_{1,0})=T_{Z^{2}}(\tilde{I})\equiv C_{Z^{2}}$. Analogous to $C_{Z}$, $C_{Z^{2}}$ is a constant matrix, i.e. it has no dependence on any $E_{i}$’s, and thus its value can be calculated beforehand. Eq. 128 implies $E_{4,2}\propto\tilde{I}=e_{4,2,1}\tilde{I}$, where $e_{4,2,1}$ is a proportionality constant. Substituting this into Eqs. 127 gives

[TABLE]

By decomposing $T_{I}=\sum^{m^{2}}_{i=1}\lambda_{i}\ket{\lambda_{i}}\bra{\lambda_{i}}$, $E_{4,1}=\sum^{m^{2}}_{i=1}e_{4,1,i}\ket{\lambda_{i}}$, $T_{Z}=\sum^{m^{2}}_{i=1}C_{Z,i}\ket{\lambda_{i}}\bra{\lambda_{i}}$, $E_{3,1}=\sum^{m^{2}}_{i=1}e_{3,1,i}\ket{\lambda_{i}}$ and $E_{2,1}=\sum^{m^{2}}_{i=1}e_{2,1,i}\ket{\lambda_{i}}$, Eq. 129 becomes

[TABLE]

Using Eq. 123, Eq. LABEL:appendix_eqn:E4_4_3 becomes

[TABLE]

where the orthogonality relation $\sum_{i^{\prime}}\left\langle\lambda_{i}|\lambda_{i^{\prime}}\right\rangle=\delta_{ii^{\prime}}$ was used in the last step. Further decomposing Eq. LABEL:appendix_eqn:E4_4_4 into a term parallel to the identity and terms perpendicular to the identity,

[TABLE]

The terms parallel to the identity in Eq. 132 are

[TABLE]

Multiplying both side of Eq. 133 by $\rho_{R}$ and taking the trace gives

[TABLE]

Substituting in Eq. 121 gives

[TABLE]

Eq. 135 is the contribution to the expectation value per site coming from polynomial degree $p=2$ i.e. the degree of $L^{2}$ with coefficient $E_{4,2}$, cf. the last term on the right hand side of Eq. 124. Comparing this result with the right hand side of Eq. 30 reveals that $e_{4,2,1}$ is $\kappa_{1}^{2}$. The terms perpendicular to the identity in Eq. 132 are not needed for the subsequent calculation of the $L^{0}$ terms. Thus there is no need to determine them.

The same procedure to obtain Eq. 135 is now applied to Eq. 126 to obtain the contribution to the expectation value from the $L^{0}$ term. Using $E_{4,2}=e_{4,2,1}\tilde{I}$, Eq. 126 becomes

[TABLE]

Using the eigenvalue decomposition above in addition to $C_{Z^{2}}\equiv T_{Z^{2}}(\tilde{I})=\sum_{i=1}^{m^{2}}C_{Z^{2},i}\ket{\lambda_{i}}$ where $C_{Z^{2},i}$ are the elements of the constant matrix $C_{Z^{2}}$, Eq. 136 becomes

[TABLE]

Using Eq. 123, Eq. LABEL:appendix_eqn:E4_3_3 becomes

[TABLE]

Separating the term parallel to the identity from terms perpendicular to it,

[TABLE]

Multiplying both sides of Eq. 139 by $\rho_{R}$ and taking the trace,

[TABLE]

Using $\ket{\lambda_{i}}\rho_{R}=0$ for the terms perpendicular to the identity and $\tilde{I}\rho_{R}=1$ for the terms parallel to the identity, Eq. 140 becomes

[TABLE]

Substituting in Eq. 135 gives

[TABLE]

The elements $e_{3,0,1}$ are to be obtained from the linear solver in Eq. 122. Eq. 142 is the contribution coming from polynomial degree $p=1$ i.e. the degree of $L^{1}$ with coefficient $E_{4,1}$, cf. Eq. 124. Comparing this result with the right hand side of Eq. 30 reveals that $e_{4,1,1}$ is $\kappa_{2}$.

VIII.2 Inset figures of the Binder cumulant

This section displays the larger version of the three inset figures corresponding to the Binder cumulant $U_{4}$ of the 1D TFI model and the 1D TKI for the purpose of clearly seeing the imperfect intersection of $U_{4}$ between the different values of $m$. Fig. 25 corresponds to the inset in the top figure of Fig. 2 corresponding to the 1D TFI model when $s=2\neq s^{*}$. At $B=1$, an imperfect intersection of $U_{4}$ between the different values of $m$ is apparent.

Figs. 26 and 27 correspond to the inset in the top and bottoms figures of Fig. 8 respectively for the 1D TKI. When $s=2\neq s^{*}$, multiple spurious crossing points occur in fig. 26. When $s=s^{*}=5.29$, the point at $J_{\perp}=2.22$ appears to be a crossing point, however upon closer inspection as shown in fig. 27, it is not - there is only one crossing point at $J_{\perp}=2.21$ as shown in the bottom figure of fig. 8.