Phase transition of four-dimensional lattice $\phi^4$ theory with tensor renormalization group

We investigate the phase transition of the four-dimensional single-component $\phi^4$ theory on the lattice using the tensor renormalization group method. We have examined the hopping parameter dependence of the bond energy and the vacuum condensation of the scalar field $\langle\phi\rangle$ at a finite quartic coupling $\lambda$ on large volumes up to $V=1024^4$ in order to detect the spontaneous breaking of the $\mathbb{Z}_2$ symmetry. Our results show that the system undergoes the weak first-order phase transition at a certain critical value of the hopping parameter. We also make a comparative study of the three-dimensional $\phi^4$ theory and find that the properties of the phase transition are consistent with the universality class of the three-dimensional Ising model.


I. INTRODUCTION
The issue of the triviality of the four-dimensional (4d) φ 4 theory has been a theoretical concern among particle physicists, because it is related to the scalar sector in the standard model [1][2][3][4][5][6][7][8][9][10][11][12][13]. The single-component φ 4 theory becomes equivalent to the Ising model in the infinite limit of the quartic coupling λ = ∞ so that numerical studies of the 4d Ising model have been performed as a nonperturbative test of the triviality, assuming the universality [14][15][16][17][18][19][20]. 1 So far, no Monte Carlo calculation has confirmed the logarithmic correction to the mean-field exponents in the scaling behavior of the specific heat, which is expected from the perturbative renormalization group analysis [26]. Moreover, a detailed Monte Carlo study has found a serious finite-volume effect due to nontrivial boundary effects in the 4d Ising model [20].
Recently, the authors have investigated the phase transition of the 4d Ising model with the higher-order tensor renormalization group (HOTRG) algorithm [27]. The tensor renormalization group (TRG) method, 2 which contains the HOTRG algorithm, has several superior features over the Monte Carlo method. (i) Since the TRG provides a deterministic numerical method, it does not have the sign problem encountered in stochastic methods, including the standard Monte Carlo simulation, as confirmed in various studies of quantum field theories [32,[35][36][37][38][39][40][41]. (ii) Its computational cost depends on the system size only logarithmically. (iii) The computational cost to simulate fermions is almost equivalent to that to bosons because the TRG can directly manipulate the Grassmann variables [32][33][34]42]. (iv) We can obtain the partition function or the path-integral itself. Thanks to the above feature (ii), we have been allowed to enlarge the lattice volume up to V = 1024 4 , which is essentially identified as the thermodynamic limit, and found finite jumps for the internal energy and the magnetization as functions of temperature in the 4d Ising model [27]. These are characteristic features of the first-order phase transition. Having shown that the 4d Ising model undergoes the weak first-order phase transition, our interest turns to the order of the phase transition in the 4d single-component φ 4 theory, which has the global Z 2 symmetry as with the Ising model. 3 In this paper, we investigate the phase transition of the 4d single-component φ 4 theory with the quartic coupling λ and the hopping parameter κ, employing the anisotropic TRG (ATRG) algorithm [30], which was proposed to reduce the computational cost of the TRG method. The ATRG has been successfully applied to analyze the 4d complex φ 4 theory at the finite density with parallel computation [41]. Our main purpose is to determine the order of the phase transition by examining the κ dependence of the bond energy and the vacuum condensation of the scalar field φ around the critical value of κ c for the fixed λ, the latter of which is an order parameter of the phase transition caused by the spontaneous Z 2 symmetry breaking. We study the model with a single choice of λ = 40, which is a finite-λ generalization of the Ising model study performed in Ref. [27], corresponding to λ = ∞. The choice of λ = 40 may also be helpful to avoid the weak coupling region affected by the Gaussian fixed point at λ = 0. For comparison, we also make the same analysis of the 3d single-component φ 4 theory at λ = 40, which is believed to belong to the universality class of the 3d Ising model. We discuss the differences between the results of the 3d and 4d cases. This paper is organized as follows. In Sec. II we explain the formulation of the lattice φ 4 theory and the ATRG algorithm. We present numerical results for the 4d and 3d cases in Sec. III and discuss the properties of the phase transition. Section IV is devoted to summary and outlook.

II. FORMULATION AND NUMERICAL ALGORITHM
We use the following popular action for the ddimensional single-component φ 4 theory on a lattice Γ: whereν is the unit vector of the ν-direction. This formulation, which is explicit about the relation to the Ising model, is equivalent to the more conventional expression The partition function is defined by using the action of Eq. (1) with the path integral measure We express the partition function as a tensor network in the similar way to Ref. [41]. The continuous variables φ n are discretized by the K-point Gauss-Hermite quadrature rule as where φ α and ω α are the α-th node and its weight. The partition function is thus discretized as where Each matrix M is approximated by the singular value decomposition (SVD) with a bond dimension D as where σ k is the k-th singular value sorted in the descending order, and U, V are the orthogonal matrices composed of the singular vectors. One finally obtains a tensor net-work representation for Z(K) as where with the shorthand notations such as i ν = i ν,n and i ν = i ν,n−ν .
In this study, we employ the parallelized d-dimensional ATRG algorithm developed in Refs. [41,47]. We keep the bond dimension D fixed throughout the ATRG procedure. For the swapping bond parts explained in Refs. [30,48], the randomized SVD is applied with the choice of p = 2D and q = 2D, where p is the oversampling parameter and q is the numbers of QR decomposition.

A. 4d case
The partition function of Eq. (12) is evaluated using the ATRG algorithm on lattices with the volume V = L 4 (L = 2 m , m ∈ N) employing the periodic boundary conditions for all the space-time directions. As explained in the previous section, there are two important algorithmic parameters. One is the number of nodes K in the Gauss-Hermite quadrature method to discretize the scalar field. The other is the bond dimension D. We check the convergence behavior of the free energy as a function of K and D by defining the following quantities: and Figure 1 shows the K dependence of δ K with D = 50 on V = 1024 4 at κ = 0.0763059 and 0.0765000, which are in the symmetric and broken symmetry phases. Note that κ = 0.0763059 is close to the transition point κ c , as we will see below. We observe that δ K decreases monotonically as a function of K and reaches the order of 10 −7 around K = 1500. This shows that the Gauss-Hermite quadrature method is not affected by whether the system is in the symmetric or broken symmetry phase. We also plot the D dependence of δ D in Fig. 2, which shows the fluctuation of free energy is suppressed as δ D ≈ 10 −5 up to D = 50. Since the double-well potential in the φ 4 theory becomes sharper for larger λ, we take a large number of K to achieve good convergence for δ K . In the following, numerical results at λ = 40 are presented for K = 2000 and D = 50 which are large enough in this study.
The phase transition point κ c is determined by following the method employed in the Ising case [27]. Suppose we have obtained a coarse-grained tensor T after the m times of coarse-graining. Defining a D × D matrix as  we calculate This quantity, introduced in Ref. [49], possibly counts the number of the largest singular value of A (m) . Therefore, it is expected that X (m) = 1 holds for the symmetric phase and X (m) = 2 for the broken symmetry phase. We may distinguish both phases by observing the plateau of X (m) after sufficient coarse-graining iterations.
In order to check the applicability of the above method to determine the value of κ c , we calculate κ c at λ = 5 and compare it with the previous results obtained by various methods including the Monte Carlo simulation [50]. Since we have found that the convergence of the free energy with respect to the bond dimension at λ = 5 becomes slightly slower than that at λ = 40, we have taken D = 55 (and K = 2000) to evaluate κ c at λ = 5. Up to D = 55, the relative error for the free energy is suppressed to O(10 −5 ). Figure 3 shows the m dependence of the value of X (m) at κ = 0.089225 and 0.089300, whose difference ∆κ = 7.5 × 10 −5 is the finest resolution across the transition point. We find X  Table III in Ref. [50]. For details on the dynamical or effective mean field theory, see Ref. [50]. For Kikuchi's method, see Ref. [51].
the critical kappa κ c = 0.0892625(375) on the 1024 4 lattice, whose error bar is provided by the resolution of κ. In Fig. 4 we find that our result is comparable to the Monte Carlo result κ c = 0.08893 (20) in Ref. [50]. Slight deviation from the Monte Carlo result may be attributed to the finite size effect: our result is obtained on the 1024 4 lattice, while the previous one is on the 32 4 lattice.
Having confirmed the validity of the method using X (m) , we determine κ c at λ = 40 with D = 50 and K = 2000. The result is κ c = 0.076305975(25) on the 1024 4 lattice, whose error bar is provided by the resolution of κ. In Fig. 5 we check the 1/λ dependence of κ c toward the Ising limit, where the result at λ = 100 is obtained in the same way as the λ = 40 case with D = 50 and K = 2000. We observe that the value of κ c seems monotonically approaching that in the Ising case. The error bars are provided by the resolution of κ but they are all within symbols.
We now turn to the investigation of the phase transi- tion with the bond energy defined by and the vacuum condensation of the scalar field φ . Both quantities are evaluated with the impure tensor method. Figure 6 plots the bond energy as a function of κ on the 1024 4 lattice. The resolution of κ becomes finer toward the transition point and the finest one is ∆κ = 5.0 × 10 −8 around the transition point. The phase transition point is consistent with κ c (gray band) determined by X (m) . Inset graph in Fig. 6 shows an emergence of a finite gap with mutual crossings of curves for different volumes, m ≥ 23, around κ c . These are characteristic features of the first-order phase transition as discussed in Ref. [52]. As the gap, we obtain by the linear extrapolation toward the transition point both from the symmetric and broken symmetry phases. In this extrapolation, we have used data points in [0.07630560, 0.07630595] for the symmetric phase and [0.0763060, 0.0763064] for the broken symmetry one. Note that we do not extrapolate ∆E b to the D → ∞ limit in this paper because a systematic study of the D dependence demands enormous computational cost and the theoretical formula for the extrapolation is not known so far. The value of ∆E b becomes smaller than the latent heat ∆E = 0.0034 (5) found in the Ising case with the HOTRG [27]. Another quantity to detect the phase transition is the vacuum condensation of the scalar field φ , which is the order parameter of spontaneous breaking of the Z 2 symmetry. We calculate φ by introducing the external fields of h = 1.0×10 −10 and 2.0×10 −10 at each κ. After taking the infinite volume limit, we extrapolate the value of φ at h = 0. Figure 7 shows the κ dependence of φ h=0 . The resolution of κ is the same as that in Fig. 6. We find that the value of κ c , where the vacuum condensation sets in, is consistent with both estimates by X (m) and the bond energy. A finite jump in φ h=0 at κ c is another indication of the first-order phase transition. We find as the value of finite jump, where we have used data points in [0.07630560, 0.07630595] for the symmetric phase and [0.0763060, 0.0763064] for the broken symmetry one, as in the case with the bond energy, to extrapolate linearly the values of φ h=0 toward the transition point. Note that this quantity is estimated as 0.037 (2) in the Ising case with the HOTRG [27].

B. 3d case
The 2d single-component lattice φ 4 theory is believed to belong to the same universality class as the 2d Ising model. The previous TRG analysis, which was carried out by two of the authors and collaborators, supports this ansatz [53]. Although the 3d case should undergo the second-order phase transition belonging to the universality class of the 3d Ising model, the direct check with the TRG method has not been performed so far. Here it must be instructive to repeat the same TRG calculation for the 3d case and compare the results between the 3d and 4d cases at λ = 40.
We first show the convergence behavior of the free energy as a function of K and D by defining the relative error in the following way: The K dependence of δ K with D = 90 on V = 4096 3 at κ = 0.112859 and 0.112920 in Fig. 8. κ = 0.112859 is near the transition point in the symmetric phase, while κ = 0.112920 is in the broken symmetric phase. We observe a monotonic decrease of δ K as a function of K, which is quite similar to the 4d case. Figure 9 shows the D dependence of δ D , where δ D reaches the order of 10 −5 up to D = 90. Notice that the achieved order of δ D is similar with the 4d case. In the following, we present the results at λ = 40 for K = 2000 and D = 90. given in the same way as Eq. (17), defining the threedimensional counterpart of Eq. (16). The value of the bond energy evaluated at κ = 0.11285900 is located within this gray band. This is due to the situation that X (m) at κ = 0.11285900 does not show any clear plateau at X (m) = 1 or 2. We observe that the bond energy on all the volumes smoothly varies as a function of κ without generating any gap. In addition, we find no mutual crossing of curves for different volumes around the phase transition point: The curve of the bond energy monotonically approaches that on the largest volume of 4096 3 . These behaviors, which are in clear contrast to the 4d case, are characteristics of the second-order phase transition as discussed in Ref. [52]. In Fig. 11, we show the κ dependence of φ h=0 , which is calculated in the same way as in the 4d case. The resolution of κ is the same as that in Fig. 10. In order to determine the transition point κ c and extract the critical exponent β, we make a fit of φ h=0 on 4096 3 lattice, which is essentially in the thermodynamic limit, employing the function of A(κ − κ c ) β over the range of κ ∈ [0.11285900, 0.11300000] in the broken symmetry phase. The fit results are A = 3.7(9), κ c = 0.112859 (6) and β = 0.32 (2). The value of β is consistent with re-cent estimates of β ≈ 0.3295 and 0.3264 for 3d Ising model with the HOTRG algorithm [29] and the Monte Carlo method [54], respectively. Numerical results for the bond energy and φ h=0 show consistency with the second-order phase transition in the universality class of the 3d Ising model.

IV. SUMMARY AND OUTLOOK
We have investigated the phase transition of the 4d single-component φ 4 theory at λ = 40 employing the bond energy and the vacuum condensation of the scalar field. Both quantities show finite jumps at the transition point on the extremely large lattice of V = 1024 4 , corresponding to the thermodynamic limit, and they indicate the weak first-order phase transition as found in the Ising limit [27]. This means that the single-component lattice φ 4 theory does not have a continuum limit. In the current ATRG calculation, the resulting latent heat ∆E b and the gap ∆ φ are smaller than those in the Ising case obtained by the HOTRG with D = 13. As a next step, it would be interesting to investigate the phase transition of the O(4)-symmetric φ 4 theory, which is more relevant to the SU(2) Higgs model.

ACKNOWLEDGMENTS
Numerical calculation for the present work was carried out with the supercomputer Fugaku provided by RIKEN (Project ID: hp200170) and also with the Oakforest-PACS (OFP) and the Cygnus computers under the Interdisciplinary Computational Science Program of Center for Computational Sciences, University of Tsukuba. This work is supported in part by Grants-in-Aid for Scientific Research from the Ministry of Education, Cul-