Phase transition of four-dimensional Ising model with higher-order tensor renormalization group

TL;DR

This paper uses the higher-order tensor renormalization group to study the four-dimensional Ising model, revealing a weak first-order phase transition through large-scale simulations.

Contribution

It applies HOTRG to a 4D Ising model at unprecedented lattice sizes, providing new insights into its phase transition nature.

Findings

Results indicate a weak first-order phase transition.

Simulations performed on lattices up to 1024^4.

Internal energy and magnetization data support the transition type.

Abstract

We apply the higher-order tensor renormalization group (HOTRG) to the four-dimensional ferromagnetic Ising model, which has been attracting interests in the context of the triviality of the scalar $\phi^4_{d=4}$ theory. We investigate the phase transition of this model with HOTRG enlarging the lattice size up to $1024^4$ with parallel computation. The results for the internal energy and the magnetization are consistent with the weak first-order phase transition.

Figures

Equations

$$Z_N=\sum _{\{\sigma =\pm 1\}}\left(\prod _{⟨ij⟩}{T}_{{\sigma }_i{\sigma }_j}\right)\left(\prod _i{V}_{{\sigma }_i}\right)$$

$${\mathcal{T}}_{i;x{x}^{\prime }y{y}^{\prime }z{z}^{\prime }t{t}^{\prime }}^{(0)}=\sum _{{\sigma }_i}{W}_{{\sigma }_{i}x}{W}_{{\sigma }_i{x}^{\prime }}{W}_{{\sigma }_{i}y}{W}_{{\sigma }_i{y}^{\prime }}{W}_{{\sigma }_{i}z}{W}_{{\sigma }_i{z}^{\prime }}{W}_{{\sigma }_{i}t}{W}_{{\sigma }_i{t}^{\prime }}{V}_{{\sigma }_i}$$

$$Z_N=\mathrm{Tr}\prod _{i=1}^N{\mathcal{T}}_i^{(0)},$$

$$Z_N\approx \mathrm{Tr}{\mathcal{T}}_{i=1}^{(n)}.$$

$$\mathrm{Tr}{\mathcal{T}}_{i=1}^{(n)}=\sum _{x,y,z,t}{\mathcal{T}}_{1;xxyyzztt}^{(n)}.$$

$${\mathcal{S}}_{i;x{x}^{\prime }y{y}^{\prime }z{z}^{\prime }t{t}^{\prime }}^{(0)}=\sum _{{\sigma }_i}{\sigma }_i{W}_{{\sigma }_{i}x}{W}_{{\sigma }_i{x}^{\prime }}{W}_{{\sigma }_{i}y}{W}_{{\sigma }_i{y}^{\prime }}{W}_{{\sigma }_{i}z}{W}_{{\sigma }_i{z}^{\prime }}{W}_{{\sigma }_{i}t}{W}_{{\sigma }_i{t}^{\prime }}{V}_{{\sigma }_i}.$$

$$⟨{\sigma }_i{\sigma }_j⟩=\mathrm{Tr}\left[{\mathcal{S}}_i^{(0)}{\mathcal{S}}_j^{(0)}\prod _{k\ne i,j}{\mathcal{T}}_k^{(0)}\right]/Z_N,$$

$${\mathcal{S}}_{1;x{x}^{\prime }y{y}^{\prime }}^{(1)}=\sum _{\alpha ,x_1,x_1^{\prime },x_2,x_2^{\prime }}{U}_{x{x}_1\otimes x_2}^{(1)}{\mathcal{S}}_{1;x_1{x}_1^{\prime }y\alpha }^{(0)}{\mathcal{S}}_{2;x_2{x}_2^{\prime }\alpha y^{\prime }}^{(0)}{U}_{x^{\prime }{x}_1^{\prime }\otimes x_2^{\prime }}^{(1)},$$

$${\mathcal{S}}_{1;x{x}^{\prime }y{y}^{\prime }}^{(n+1)}=\sum _{\alpha ,x_1,x_1^{\prime },x_2,x_2^{\prime }}{U}_{x{x}_1\otimes x_2}^{(n+1)}{\mathcal{S}}_{1;x_1{x}_1^{\prime }y\alpha }^{(n)}{\mathcal{T}}_{2;x_2{x}_2^{\prime }\alpha y^{\prime }}^{(n)}{U}_{x^{\prime }{x}_1^{\prime }\otimes x_2^{\prime }}^{(n+1)}.$$

$$⟨{\sigma }_i{\sigma }_j⟩\approx \frac{\mathrm{Tr}{\mathcal{S}}_1^{(n)}}{\mathrm{Tr}{\mathcal{T}}_1^{(n)}}.$$

$$⟨{\sigma }_i⟩=\mathrm{Tr}\left[{\mathcal{S}}_i^{(0)}\prod _{k\ne i}{\mathcal{T}}_k^{(0)}\right]/Z_N,$$

$$⟨{\sigma }_i⟩\approx \frac{\mathrm{Tr}{\mathcal{S}}_1^{(n)}}{\mathrm{Tr}{\mathcal{T}}_1^{(n)}}.$$

$$\delta f=\left|\frac{\ln Z_N(D_{\mathrm{cut}})-\ln Z_N(D_{\mathrm{cut}}=13)}{\ln Z_N(D_{\mathrm{cut}}=13)}\right|.$$

$$A_{t{t}^{\prime }}^{(n)}=\sum _{x,y,z}{\mathcal{T}}_{1;xxyyzzt{t}^{\prime }}^{(n)},$$

$$X^{(n)}=\frac{{\left(\mathrm{Tr}{A}^{(n)}\right)}^2}{\mathrm{Tr}{\left(A^{(n)}\right)}^2},$$

$$\Delta E(D_{\mathrm{cut}}=13)=0.0034(5),$$

$$\Delta m(D_{\mathrm{cut}}=13)=0.037(2).$$

Full text

††preprint: UTHEP-723, UTCCS-P-115

1 Introduction

It is well known that the critical behavior of the Ising model on the higher-dimensional hypercubic lattice is well explained with the mean-field theory. In dimensions larger than four, the effect of the background fluctuations becomes negligible and the model in the critical region exactly obeys the mean-field exponents Aizenman:1981zz ; Aizenman:1982ze . At the upper critical dimension, however, multiplicative logarithmic corrections are added to the leading scaling behavior of the mean-field theory. Some of these corrections were derived by the perturbative calculation with the renormalization group method PhysRevB.7.248 . Since the Ising model is specified by the infinite coupling limit of the single-component scalar $\phi^{4}_{4}$ theory, the model in four dimensions has been attracting the interest of particle physicists for a long time in the context of the triviality of the $O(4)$ scalar $\phi^{4}_{4}$ theory, which corresponds to the scalar sector of the standard model describing the generation of gauge boson and fermion mass through the Higgs mechanism Wilson:1973jj ; Luscher:1987ay ; Luscher:1987ek ; Luscher:1988uq ; Huang:1988hu ; Jansen:1988cw ; Frick:1989gw . There is also a recent study to discuss the triviality of the $O(N)~{}\phi^{4}_{4}$ theory with the higher-loop beta function Shrock:2014zca ; Shrock:2016hqn ; Shrock:2017zuk .

A numerical study of the Ising model on a hypercubic lattice serves as a non-perturbative test of the triviality Kenna:1992np ; Kenna:2004cm ; if the leading scaling behavior is the mean-field type and it is modified only by the multiplicative logarithmic factor, one obtains a supporting evidence for the triviality. In fact, the numerical investigation based on the Monte Carlo simulation has successfully caught the mean-field exponents PhysRevB.22.4481 ; SanchezVelasco:1987ah ; Bittner:2002pk ; PhysRevE.80.031104 ; Lundow:2010en , but there remains some controversy over the appearance of the logarithmic corrections PhysRevE.80.031104 ; Lundow:2010en ; Cea:2005ad ; Stevenson:2005yn ; Balog:2006fs . Actually no Monte Carlo study has confirmed the logarithmic correction in the scaling behavior of the specific heat, which is $(\ln|t|)^{1/3}$ with $t$ the reduced temperature expected from the perturbative renormalization group analysis. This is mainly because the cubic root of logarithmic divergence is too weak to detect by the finite size scaling analysis, or the specific heat may be actually bounded PhysRevE.80.031104 . Indeed, the finite volume effect of the four-dimensional Ising model had been investigated from various viewpoints Jansen:1988cw . A detailed Monte Carlo study has found serious finite-volume effect due to non-trivial boundary effects in the four-dimensional Ising model Lundow:2010en . From the viewpoint of numerical calculation, it could be possible that there remain some unrevealing aspects in the phase transition of this model and it should be worth trying different approaches other than the Monte Carlo method.

For this purpose we employ the tensor network scheme to investigate the four-dimensional classical Ising model. This scheme has various types of numerical algorithms Orus:2018dya , which can be divided into two streams: Hamiltonian approach and Lagrangian one. The latter enables us to evaluate the partition functions directly via tensor network representation. A typical algorithm is the tensor renormalization group (TRG) Levin:2006jai , which was originally proposed by Levin and Nave for the two-dimensional Ising model. The TRG method has been successfully applied to the two-dimensional field theories with the path-integral formulation in the particle physics Shimizu:2012zza ; Shimizu:2012wfa ; Liu:2013nsa ; Denbleyker:2013bea ; Shimizu:2014uva ; Shimizu:2014fsa ; Unmuth-Yockey:2014afa ; Takeda:2014vwa ; Kawauchi:2016dcg ; Meurice:2016mkb ; Shimizu:2017onf ; Kadoh:2018hqq ; Sakai:2018xkx ; Kadoh:2018tis . The higher-order TRG (HOTRG) PhysRevB.86.045139 is an improvement of TRG with the extension to higher dimensions. One of attractive features in TRG and HOTRG is that we are allowed to directly study the thermodynamic properties; we can systematically increase the system size by repeating the coarse-graining steps in the algorithms. Although earlier studies with HOTRG are restricted to two- and three-dimensional systems Qin_2013 ; Yu:2013sbi ; Wang_2014 ; Kawauchi:2015heu ; Genzor:2015pua ; PhysRevB.93.125115 ; Kawauchi:2016xng ; Sakai:2016tzv ; Sakai:2017jwp ; Chen:2017ums ; Yoshimura:2017jpk ; PhysRevE.98.062114 ; Chen:2018nzs ; Kuramashi:2018mmi ; NISHINO2019 , including the three-dimensional classical Ising model PhysRevB.86.045139 , the algorithm itself is readily extended to a four-dimensional lattice. In this paper, we employ the HOTRG method to investigate the phase transition of the classical Ising model on the four-dimensional hypercube. The accuracy of HOTRG is controlled by the bond dimension $D_{\mathrm{cut}}$, which is varied up to 14 in this study. In order to investigate the phase transition of the model, we measure the internal energy and the magnetization through the evaluation of the tensor network with some impurity located at the center of hypercube.

This paper is organized as follows. In Sec. 2 we briefly review the HOTRG method and explain an approach to evaluate the internal energy and the magnetization of the Ising model. We present numerical results in Sec. 3 and discuss the properties of the phase transition. Sec. 4 is devoted to summary and outlook.

2 HOTRG w/ and w/o impurity

The partition function of the four-dimensional ferromagnetic Ising model is given by

[TABLE]

with $T_{\sigma_{i}\sigma_{j}}=\mathrm{e}^{\beta\sigma_{i}\sigma_{j}}$, $V_{\sigma_{i}}=\mathrm{e}^{\beta h\sigma_{i}}$ , where $\sigma_{i}$ is the two-state classical spin variable on the lattice site $i$, $\braket{ij}$ specifies the sum over all the nearest-neighboring spin pairs, $\beta$ is the inverse temperature $1/T$ and $h$ is the external magnetic field. The subscript $N$ is the size of a system. Based on the eigenvalue decomposition $T=U\Lambda U^{\mathrm{T}}$, one defines the eight-rank local tensor locating on each lattice site as

[TABLE]

where $W=U\sqrt{\Lambda}$. The indices of these tensors are called bond indices. Now, we obtain the tensor network representation of Eq. (1) as

[TABLE]

where we assume the periodic boundary condition and the right-hand side means all the bond indices are contracted so as to restore the model defined on the four-dimensional hypercube. One way to evaluate Eq. (3) is HOTRG with the use of the higher-order singular value decomposition (HOSVD) PhysRevB.86.045139 . In the HOTRG procedure, nearest two local tensors along the $x$-, $y$-, $z$- and $t$-directions are mapped to the coarse-grained one sequencially. Hence the lattice size is reduced by a factor of 2 after each step of coarse-graining. After repeating $n$ steps of coarse-graining, one obtains the partition function with the system size of $N=2^{n}$; that is,

[TABLE]

The right-hand side is again the sum over all the bond indices so as to restore the structure of the four-dimensional lattice model with the periodic boundary condition and this is easily done by defining the trace of the coarse-grained tensor as

[TABLE]

There are two ways to evaluate the expectation values such as the internal energy and the magnetization. One is the numerical differentiation with respect to $\beta$ and $h$. The other is the direct evaluation of the expectation value using the corresponding tensor network representation. For instance, we can obtain the internal energy through the evaluation of the nearest-neighbor local energy term $\Braket{\sigma_{i}\sigma_{j}}$ with the HOTRG method as follows. We first define the additional local tensor as

[TABLE]

With the use of this local tensor, the tensor network representation for the local energy is given by

[TABLE]

where $\mathcal{S}_{i,j}^{(0)}$ represent the tensors on the lattice sites $i$ and $j$, respectively, and $\mathcal{T}_{k}^{(0)}$ is for the rest of $N-2$ sites. Since the numerator looks as if it contains two impurities, we call $\mathcal{S}_{i}^{(0)}$ impure tensor and $\mathcal{T}_{k}^{(0)}$ pure tensor. The denominator is evaluated by the plain HOTRG method. In order to coarse-grain the impure tensor network of Eq. (7), we assume the local energy term is fixed at the center of lattice ($i=1$, $j=2=1+{\hat{y}}$ with ${\hat{y}}$ the unit vector in $y$-direction) during the HOTRG calculation. At the first step, we define the coarse-grained impure tensor $\mathcal{S}_{1}^{(1)}$ by contracting two initial impurities. For simplicity we give the corresponding expression in the two-dimensional case;

[TABLE]

where $U^{(1)}$ is a block-spin transformation determined within the original algorithm of HOTRG PhysRevB.86.045139 . In the following steps, $\mathcal{S}^{(n+1)}_{1}$ is defined by the combination of $\mathcal{S}_{1}^{(n)}$ and $\mathcal{T}_{2}^{(n)}$. We again show the corresponding formula in two-dimensional case for simplicity;

[TABLE]

Finally, the local energy is approximately given by

[TABLE]

The meaning of the trace is the same as in Eq. (5). Since the original model has the translational invariance, $\Braket{\sigma_{i}\sigma_{j}}\times d$, where $d$ is the dimensionality, should give the absolute value of internal energy.

The one-point function $\Braket{\sigma_{i}}$ to measure the magnetization is also evaluated in the same way. In this case, we are allowed to apply the four-dimensional counterpart of Eq. (9) from the first coarse-graining step because the initial expression of $\Braket{\sigma_{i}}$ has the form

[TABLE]

where $\mathcal{S}_{i}^{(0)}$ locates only on the lattice site $i$ ($i=1$) and $k$ runs the rest of $N-1$ sites. After sufficient iterations, $\Braket{\sigma_{i}}$ is evaluated by

[TABLE]

Thanks to the translational invariance, $\Braket{\sigma_{i}}$ directly corresponds to the spatial average of the Ising spin. Note that Eqs. (10) and  (12) have the same expression, but they are evaluated by different coarse-graining procedures.

In Ref. PhysRevB.86.045139 , computational costs and memory space requirements in two- and three-dimensional HOTRG are given. Computational costs are $\mathcal{O}(D_{\mathrm{cut}}^{7})$ and $\mathcal{O}(D_{\mathrm{cut}}^{11})$ and memory space requirements are $\mathcal{O}(D_{\mathrm{cut}}^{4})$ and $\mathcal{O}(D_{\mathrm{cut}}^{6})$, respectively. In straightforward expansion of the HOTRG algorithm in Ref. PhysRevB.86.045139 to four dimensions, the computational cost is $\mathcal{O}(D_{\mathrm{cut}}^{15})$ and the memory space requirement is $\mathcal{O}(D_{\mathrm{cut}}^{8})$. In our implementation, the computational cost in each process is $\mathcal{O}(D_{\mathrm{cut}}^{13})$ and the memory space requirement in each process is $\mathcal{O}(D_{\mathrm{cut}}^{7})$. Reduction of the order of computational cost is achieved by using $D_{\mathrm{cut}}^{2}$ processes. This implementation is basically based on an idea which will be shown in Ref. Yamashita: .We have carried out a detailed measurement of the internal energy and the magnetization with $D_{\mathrm{cut}}=13$ employing the fine resolution of the temperature $\Delta T=6.25\times 10^{-6}$ around the transition temperature. We have repeated the calculation with $D_{\mathrm{cut}}=14$ to confirm the qualitative features obtained with $D_{\mathrm{cut}}=13$. But, in this case the temperature resolution remains coarser as $\Delta T=3.0\times 10^{-5}$ due to the computational cost. In the following we focus on the results with $D_{\mathrm{cut}}=13$.

As found in Sec. 3 all the measured physical quantities seem to lose the volume dependence beyond $n\approx 30$ so that the lattice size of $N=2^{40}=1024^{4}$ is large enough to be taken as the thermodynamic limit.

3 Numerical results

We first evaluate the free energy with the plain HOTRG method. The convergence behavior is investigated by defining the following quantity;

[TABLE]

Figure 1 shows a typical convergence behavior of $\ln Z_{N}$ in the vicinity of the transition temperature. We observe that $\delta f$ decreases monotonically as a function of $D_{\mathrm{cut}}$.

We now turn to the determination of the transition temperature. Let us assume that one has just obtained the coarse-grained tensor $\mathcal{T}^{(n)}_{1;xx^{\prime}yy^{\prime}zz^{\prime}tt^{\prime}}$ whose coarse-graining direction was $t$-direction. Choosing a $D_{\mathrm{cut}}\times D_{\mathrm{cut}}$ matrix as

[TABLE]

we calculate the following quantity;

[TABLE]

which counts the number of the largest singular value of $A^{(n)}$. This is an indicator of the symmetry-breaking PhysRevB.80.155131 . We calculate $X^{(n)}$ iteratively until it converges. A typical convergence behavior of $X^{(n)}$ is shown in Fig. 2. Notice that we sequentially redefine $A^{(n)}$ corresponding to the direction of coarse-graining in the practical calculation. Figure 3 shows the transition temperature $T_{\mathrm{c}}$ as a function of $D_{\mathrm{cut}}$. The error bars, provided by the temperature resolution, are all smaller than the corresponding symbols. Since $T_{\mathrm{c}}(D_{\mathrm{cut}})$ is estimated by $X^{(n)}$ with sufficiently large $n$, typically beyond $n=30$, there remains little finite-volume effect. In this work, we have obtained $T_{\mathrm{c}}(D_{\mathrm{cut}}=13)=6.650365(5)$ on $1024^{4}$ lattice. The recent Monte Carlo study PhysRevE.80.031104 obtained $\beta_{\mathrm{c}}=0.1496947(5)$ corresponding to $T_{\mathrm{c}}=6.68026(2)$, which shows a slight deviation from our result with HOTRG up to $D_{\mathrm{cut}}=13$. Note that the value of $T_{\mathrm{c}}$ in Ref. PhysRevE.80.031104 was obtained by the infinite-volume extrapolation using the results on relatively small lattices with $L^{4}\leq 80^{4}$.

Let us move on to the evaluation of the internal energy, which can be obtained by numerical differentiation or the coarse-graining of the impure network of Eq. (7). We have compared both methods varying the temperature resolution and found that the latter successfully keeps the numerical accuracy as the resolution becomes finer. In the following, we show the results with the impure tensor method. Figure 4 traces the volume dependence of the internal energy with $D_{\mathrm{cut}}=13$. The converging behavior toward the thermodynamic limit is clearly observed. Since the system size $N$ is given by $2^{n}$, a hypercubic structure is restored in the condition of $n~{}\mathrm{mod}4=0$. Figure 5 shows the internal energy as a function of temperature for various lattice sizes with $D_{\mathrm{cut}}=13$. In the case of $n\geq 24$ ($L\geq 64$), a finite jump emerges with mutual crossings of curves between different volumes around the transition temperature. These are characteristic features of the first-order phase transition as discussed in Ref. Fukugita1990 . The similar volume dependence and a finite jump at $L\geq 64$ have been also confirmed in case of $D_{\mathrm{cut}}=14$. The numerical value of the finite jump $\Delta E(D_{\mathrm{cut}}=13)$ in the infinite-volume limit is

[TABLE]

which is obtained by the linear extrapolation toward the transition temperature both from the low and high temperature regions. The resolution of the temperature at the boundary between the two phases is $\Delta T=6.25\times 10^{-6}$.

We also investigate the spontaneous magnetization, which is an order parameter to detect the symmetry-breaking phase. Figure 6 shows a typical volume dependence of magnetization toward the thermodynamic limit. We have evaluated $\Braket{\sigma_{i}}$ with $h=1.0\times 10^{-9}$ and $2.0\times 10^{-9}$ at each temperature and coarse-graining step. After taking the infinite-volume limit we extrapolate the value of $\Braket{\sigma_{i}}$ toward the $h\to 0$ limit. Figure 7 shows the resulting spontaneous magnetization as a function of temperature. The transition temperature is consistent with both estimates by $X^{(n)}$ and the internal energy. We have observed a finite jump in the magnetization, whose numerical value is obtained by the linear extrapolation toward the transition temperature both from the low and high temperature regions;

[TABLE]

The resolution of the temperature at the boundary between the two phases is again $\Delta T=6.25\times 10^{-6}$. Note that we have tried several choices of the external field other than $h=\mathcal{O}(10^{-9})$ and confirmed that the behavior of the magnetization is robust against the change of the magnitude of $h$.

4 Summary and outlook

We have analized the phase transition of the four-dimensional ferromagnetic Ising model employing HOTRG on $L^{4}\leq 1024^{4}$ lattices. The transition temperature is successfully determined by measuring the degeneracy of the largest singular value of the pure tensor. We have also investigated the temperature dependence of the internal energy and magnetization with the impure tensor method. We have found a finite jump for the internal energy together with mutual crossings of curves between different volumes around the transition temperature. A finite jump is also observed in the magnetization. These are characteristic features of the first-order phase transition. The numerical results obtained by the impure tensor method are consistent with the weak first-order phase transition. The resulting estimate for the transition temperature in the thermodynamic limit shows a slight deviation from the recent Monte Carlo prediction PhysRevE.80.031104 obtained from the infinite-volume extrapolation of the data on relatively small lattices up to $80^{4}$.

In future investigation, the HOTRG calculation with $D_{\mathrm{cut}}>14$ should allow us to achieve a direct and essential improvement of this study. Our impure tensor method can be also improved by considering all the patterns of coarse-graining for the network including some impurities Morita:2018tpw . Another possible approach is to develop the best optimization of the Frobenius norm of impure tensor, which would be a realistic way to improve our impure tensor algorithm from the viewpoint of the computational cost of HOTRG.

Acknowledgements.

Numerical calculation for the present work was carried out with the COMA (PACS-IX) computer under the “Interdisciplinary Computational Science Program" of Center for Computational Sciences, University of Tsukuba and with the Oakforest-PACS system of Joint Center for Advanced High Performance Computing. This work is supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) as “Exploratory Challenge on Post-K computer (Frontiers of Basic Science: Challenging the Limits)”.