Tensor network analysis of critical coupling in two dimensional ϕ4 theory

We make a detailed analysis of the spontaneous Z2-symmetry breaking in the two dimensional real ϕ4 theory with the tensor renormalization group approach, which allows us to take the thermodynamic limit easily and determine the physical observables without statistical uncertainties. We determine the critical coupling in the continuum limit employing the tensor network formulation for scalar field theories proposed in our previous paper. We obtain [λ/μc2]cont. = 10.913(56) with the quartic coupling λ and the renormalized critical mass μc. The result is compared with previous results obtained by different approaches.


Introduction
The two dimensional real φ 4 theory has been extensively studied as the simplest model with a spontaneous symmetry breaking which plays important roles in particle physics. Although the continuous symmetry as in the complex φ 4 theory cannot be broken in two dimensions [1], the real φ 4 theory has a discrete Z 2 -symmetry that can be spontaneously broken [2]. The critical coupling, which separates two phases, has been studied by various approaches with and without lattice simulations . While many Monte Carlo studies have been carried out, there are several studies using deterministic schemes such as matrix product state (MPS), Hamiltonian truncation methods, and Borel resummation. In both directions, the thermodynamic limit is important to seriously study the critical phenomena.
The tensor renormalization group (TRG) [27] enables us to take the thermodynamic limit easily: the computational cost of the TRG depends on the space-time volume logarithmically while that of the Monte Carlo methods is proportional to the volume. Furthermore, the TRG method does not contain any stochastic procedure to evaluate physical quantities. Then, using the TRG, we can observe precisely the symmetry breaking effects with larger volumes avoiding statistical uncertainties.
Since the TRG was originally proposed for the two dimensional Ising model, several improved algorithms [28][29][30] and extensions to fermion systems [31,32] have been developed. Those techniques are introduced into the path-integral formulation of the field theory JHEP05(2019)184 and many studies have been carried out [33][34][35][36][37][38][39][40][41][42][43][44][45]. 1 It is, however, not straightforward to apply the TRG to lattice scalar theories because the indices of tensor are essentially the field variables so that their dimension is infinite unlike the Ising case. 2 An attempt to construct a tensor network (TN) in scalar field theory was made in ref. [47], where an orthogonal function expansion was employed to define a finite dimensional tensor. It, however, suffers from unavoidable loss of significance through its numerical procedure. To overcome this difficulty, the authors have recently proposed another formulation in which the Gaussian quadrature formula is used to discretize the field variables [48]. Although such a simple discretization of continuous variables has already been attempted in creating transfer matrices [49][50][51][52][53][54], we use it to represent the partition function as a tensor network in two dimensional system in ref. [48]. The effectiveness of our formulation has been confirmed only in a non-interacting system. In the current paper, with studying the critical phenomena, we also aim to establish our TN formulation for an interacting case.
In this paper we determine the critical coupling of the two dimensional lattice φ 4 theory employing the TRG method with our TN formulation. The critical coupling is determined from the power law behavior of the susceptibility for the expectation value of the scalar field near the critical point. We extrapolate the critical coupling to the continuum limit and compare it to other recent results obtained by different approaches [16,17,20,26]. We also make some analyses to estimate the magnitude of systematic error from the discretization of the field variables.
This paper is organized as follows. We present the two dimensional continuum and lattice φ 4 theories to fix our notation in section 2. A TN formulation of the expectation value of the scalar field is given in section 3. In section 4, numerical results are presented. Some details of our method and systematic errors are discussed in section 4.1. The results of susceptibility are then presented in section 4.2, and the critical couplings at several lattice spacings are listed in table 1 of section 4.3. We show the final result of the critical coupling in the continuum limit in section 4.4. Section 5 is devoted to summary and discussion. In appendix A the coarse-graining algorithm of TN with an impurity tensor is explained. The systematic error originated from the Gaussian quadrature is discussed in appendix B.

Two dimensional φ 4 theory
The Euclidean continuum action of the two dimensional real φ 4 theory is defined as with a real scalar field φ (x) ∈ R, the bare mass µ 0 , and the quartic coupling constant λ > 0. This theory is super renormalizable. The mass renormalization is required only 1 Separately from the path-integral formulation employed in this paper, Hamiltonian based approaches also show a remarkable developments; for the detail see a recent review given in ref. [46]. 2 In lattice gauge theories, one can use the character expansion with a truncation to construct a finite dimensional tensor.

JHEP05(2019)184
at one-loop level while the renormalization of the coupling constant is not necessary to all orders of perturbation theory. The suffix zero is used to indicate that the mass is a bare parameter and to be renormalized. This model has the Z 2 -symmetry (φ → −φ) classically, but it may be spontaneously broken at a quantum level. The expectation value of the scalar field φ (x) , which does not actually depend on the coordinate x from the translational invariance, is an order parameter of the symmetry breaking: φ (x) = 0 for the symmetric phase and φ (x) = 0 for the symmetry broken phase.
The dimensionful coupling λ gives a typical scale of this model while a dimensionless coupling characterizes its physical behavior, where µ is the renormalized mass. At the tree level, two phases and the sign of κ 0 = λ/µ 2 0 are in one-to-one correspondence, and the critical point is κ 0,c = 0. The negative one corresponds to the broken phase. However, the critical point can change with the quantum corrections. The lattice simulations play a crucial role to determine the renormalized critical point κ c beyond the perturbation theory.
Let us define the lattice theory on a square lattice where a positive integer L is the lattice size. The lattice spacing a is assumed to be a = 1, and we suppress it unless necessary. We label the lattice site (n 1 , n 2 ) simply as n. The lattice scalar field φ n is defined on the site n and satisfies the periodic boundary conditions, whereρ is the unit vector of the ρ-direction. The lattice action is given by Note that µ 0 and λ are dimensionless quantities written as a 2 µ 2 0 and a 2 λ, respectively. The mass renormalization can be calculated by the lattice perturbation theory. We employ the same renormalization condition as in refs. [2,16]. The renormalized mass µ is defined as which is the one-loop self energy. We use eq. (2.6) to determine the renormalized critical coupling in later section.

JHEP05(2019)184
3 Tensor network formulation for two dimensional lattice φ 4 theory The expectation value of the scalar field can be expressed as a TN form. We present the detailed procedure following our previous work [48]. We begin with the partition function The action S h is given by where a constant external field h is introduced to investigate the phase transition taking h → 0 after L → ∞. The Boltzmann weight can be expressed as a product of local factors, The function f (φ 1 , φ 2 ) is a symmetric matrix with the continuous indices φ 1 , φ 2 ∈ R, and the continuity makes numerical treatments hard.
To define a finite dimensional tensor, we use the Gauss-Hermite quadrature formula, in which the integral of a target function g (x) with the weight function e −x 2 is approximated by a discrete sum as where y α (α = 1, 2, . . . , K) is the α-th root of the K-th Hermite polynomial H K (y) and w α is the associated weight. 3 K is the order of approximation. If g (x) is a polynomial function with the degree of 2K − 1 or less, the Gauss-Hermite quadrature is known to reproduce the exact result. For general functions, the convergence of the Gauss-Hermite quadrature is not obvious, and in this paper we numerically check it for our particular case. By replacing each integral in eq. (3.1) by eq. (3.5) and using eq. (3.3), we obtain a discretized version of the partition function The n-th Hermite polynomial is defined by yα and weight wα of the n-point Gauss-Hermite quadrature are given by Hn (yα) = 0 and wα = 2 n−1 n! √ π/(n 2 Hn−1 (yα) 2 ). See ref. [55] for the detail.

JHEP05(2019)184
, it can be decomposed by the singular value decomposition (SVD) as 4 where σ i are the singular values sorted as σ 1 ≥ σ 2 ≥ · · · ≥ σ K ≥ 0 and U, V are unitary matrices. The range of the summation in eq. (3.7) will be truncated to reduce the computational cost later. To this end, we used the SVD here. Finally we have . Now let us turn to the expectation value of φ, Since Z is already expressed as a uniform TN as shown in eq. (3.8), the remaining task is to represent the numerator Z 1 as a TN. The insertion of the operator φ n in eq. (3.11) yields where an extra y α is inserted compared to eq. (3.9). The modified tensorT . is referred to as impurity tensor. Repeating the same procedures as for the partition function, we obtain the following TN form: in which one of T s in eq. (3.9) is replaced by the impurity tensorT . The contraction of tensor indices in eq. (3.8) is exactly taken from 1 to K. In actual computations, we initially truncate the bond dimension from K to D on account of computational complexity; namely we redefine {x,t} as n∈Γ L D xn=1 D tn=1 , and then the JHEP05(2019)184 partition function in eq. (3.8) is initially approximated. In TRG, renormalized tensors have a truncated bond dimension D as a tunable parameter. In this study we set D = D. 5 Now we have two parameters to control the accuracy of the approximations, the degree of the Gauss-Hermite quadrature K and the size of tensors D, and we will check the stability of numerical results with respect to D and K in appendix B.

Numerical results
In this section the numerical results are shown step by step. First, the technical details related to the Gauss-Hermite quadrature and the TRG are summarized in section 4.1. In section 4.2 the susceptibility of φ in the thermodynamic limit is evaluated. Using the susceptibility we determine the critical parameters in section 4.3, and in section 4.4 we take the limit λ → 0 to obtain the critical coupling in the continuum limit.

Methods
The expectation value φ is obtained by the ratio of Z (K) and Z 1 (K) shown as eq. (3.10), which are individually evaluated by the TRG method with the truncated bond dimension D ∈ [16,64]. One can use the standard technique for computing Z (K) in eq. (3.8) and needs a little ingenuity for Z 1 (K) in eq. (3.13). The coarse-graining procedure of the tensor network with an impurity tensor is shown in refs. [37,47,56,57], and the detailed technique is also described in appendix A.
We remark that the systematic errors can be estimated by investigating the K-and D-dependences of the numerical results. There are three sources of systematic errors in the present case: (a) the discretization of the scalar fields with the K-point Gauss-Hermite quadrature, (b) the truncation of the bond dimension of the initial tensor discussed in the last paragraph in section 3, and (c) the truncation in coarse-graining steps of TRG. The effect of (a) is appeared as K-dependence of the results while the errors from (b) and (c) are observed through the D-dependence of the results.
In figure 1 we plot the singular values of the SVD in eq. (3.7). The values of µ 0 are chosen near the criticality for a given λ. A clear hierarchy of the singular values assures that the modification of contraction range that is explained at the end of the last section effectively work. We fix K = 256 that is large enough so that the systematic errors associated with the choice of K are to be smaller than those of D as discussed in appendix B. In section 4.3 we investigate the D-dependence of the critical coupling λ/µ 2 c and estimate the systematic errors by measuring small fluctuations which occur when D changes.

Susceptibility
The critical point is given by the peak position of the susceptibility for φ defined by  .7), which is used to define the initial tensor in eq. (3.9).
where φ h,L denotes the expectation value measured for a given constant external field h and lattice size L. 6 Equation (4.1) is reduced to χ = lim h→0 lim L→∞ φ h,L /h since φ 0,L = 0. It is rather complicated procedure to take the double limits with respect to h and L numerically. We first find a constant behavior of φ h,L /h as L increases, which effectively represents φ h,∞ /h. Then, χ is obtained as a converged value of φ h,∞ taking h → 0.
In figure 2 we show φ h,L /h as a function of L for several values of h. One can see that it becomes a constant for larger volumes L ≥ 10 6 , which effectively corresponds to the thermodynamic limit within the systematic errors. 7 φ h,∞ /h is thus obtained for several small values of the external field h. As a representative case we choose D = 32, λ = 0.05 and, µ 2 0 = µ 2 0,rep. with in the symmetric phase very close to the critical point: 1 − µ 2 0,rep. /µ 2 0,c ≈ 10 −5 . Figure 3 shows the h-dependence of φ h,∞ /h with the same parameter set as in figure 2. It seems to converge to a constant for sufficiently small h ≤ 10 −11 . Figure 4 shows a zoomed version of figure 3, and there the values of φ h,∞ /h lie on a quadratic function. This is because φ h,∞ is an odd function of h as a consequence of the Z 2 6 One can evaluate the susceptibility also by differentiating the free energy twice with respect to h numerically, but it may suffer from a loss of significant digits. Since φ is obtained as a direct output of the TRG, we use eq. (4.1) to avoid such a problem. symmetry and φ h,∞ /h = c 0 + c 2 h 2 + · · · for h 1 with coefficients c i (i = 0, 2, . . .). Then we simply fit them using a quadratic function to determine χ, and its error is defined as the difference between the obtained susceptibility and φ h,∞ /h at h = 10 −11 . Similar behaviors are observed for other parameter sets employed in this work (see table 1). We obtain the susceptibility for all the parameter sets taking the same procedure as here.

Critical mass for each λ
The scaling behavior of χ near the critical point in the symmetric phase is given by with the critical bare mass µ 0,c , the critical exponent γ, and a constant A. This formula can be used to determine µ 0,c by fitting the numerical result of χ. It is useful to express eq. (4.3) as The φ 4 theory is expected to belong to the universality class of the two dimensional Ising model, whose critical exponent is exactly known to be γ Ising = 7/4 = 1.75 [58]. Equation (4.4) could make it easier to check the consistency between γ and γ Ising because χ −1/γ Ising is linear in µ 2 0 if γ = γ Ising . Figure 5 shows χ −1/γ Ising as a function of µ 2 0 at λ = 0.05 and D = 32. It is clear that the results are on a straight line which implies γ = γ Ising . We determine the values of µ 2 0,c by fitting the susceptibility using the expression of eq. (4.4) with µ 2 0,c and A as free parameters and fixed γ = γ Ising .
The bare critical mass µ 2 0,c for another value of λ is obtained by repeating the procedures above. 9 Then the renormalized critical mass µ 2 c is determined by solving the self-consistent equation (2.6) using µ 2 0,c and λ as inputs. 9 In this context the critical mass is regarded as a function of λ. To emphasize this we do not hide the argument λ in eq. (4.5). Note that the critical mass also depends on D as shown in table 1, but we extract the D-independent ones to take the continuum limit as discussed in the last of this section.  Table 1 summarizes the results of µ 2 0,c obtained for 0.005 ≤ λ ≤ 0.1 and 16 ≤ D ≤ 64. We also present the dimensionless critical coupling λ/µ 2 c . Since χ 2 /d.o.f. ≤ 1 for all λ and D, the linear fittings with fixing γ = γ Ising are reasonable. The critical bare mass µ 2 0,c monotonically decreases as λ approaches zero while it does not show any smooth dependence on D at each λ.
In figure 6 we plot the D-dependence of λ/µ 2 c for λ = 0.05. In the large D region λ/µ 2 c converges to a value while oscillating. This kind of behavior is also observed in the TRG computations of other models such as the Ising and the Potts models [59]. As also shown in figure 6 we determine the central value of λ/µ 2 c as an average between the maximum and minimum values in the oscillating region. The error is estimated using the half width of the oscillation. Although the critical coupling for each D has a rather small systematic error as given in table 1, its magnitude is negligible compared to the fluctuation of the critical coupling with the variation of D. Thus, we present the results of the renormalized critical coupling λ/µ 2 c only with the systematic errors associated with D.

Continuum limit of the critical coupling
We finally obtain the value of the critical coupling in the continuum limit as . (4.5) Note that the continuum limit is understood as a 2 λ → 0. Figure 7 shows the λ-dependence of the critical coupling with the systematic errors. Let us take the continuum limit (λ → 0) of eq.   Table 1. Bare critical mass µ 2 0,c , renormalized critical coupling λ/µ 2 c , and χ 2 /d.o.f. of the linear fitting described in the main body of the text for given λ and D. Errors are originated from that of the susceptibility as explained in section 4.2.   with the systematic error from the D-dependence in the TRG. Figure 8 shows a comparison among our result and the recent Monte Carlo results in refs. [16,17,20,26] for small λ. 10 We have reached the smallest lattice spacing: λ = 0.005. The error bar is relatively large compared to the latest Monte Carlo result around λ ≈ 0.01. However, such a direct comparison about the efficiency cannot be made because the strategy used in this paper is rather different from that of the worm algorithm. References [20,26] employ the condition L/ξ = const. with the lattice extent L and the correlation length ξ to search the critical coupling. This condition implies that ξ grows linearly with L, so that we arrive at the critical point when L/a → ∞ with the lattice spacing a. But what they have done is fixing the maximum value of L/a to 384 or 512 for all the calculation at different lattice spacings. In this case a longer extrapolation to the critical coupling is required towards the continuum limit, and the finite size corrections could contaminate the results more severely. Actually, figure 8 shows that the results of ref. [20] agree well with thoese of ref. [16] at the coarser lattice spacing while they become systematically deviated towards the continuum limit. Although further studies are thus needed for a direct comparison, the final results seem to be roughly consistent with each other.    [20], and Bronzin et al. [26]) and this work. At λ = 0, data points are horizontally shifted to ensure the visibility. Note that the results by Wozar and Wipf cannot be compared at non-zero λ since they used the SLAC derivative for scalar bosons, but, in the continuum limit (λ = 0), their results are consistent with naively discretized ones within errors.

JHEP05(2019)184 5 Summary and outlook
We have studied the two dimensional lattice φ 4 theory using the TRG with our new formulation of making a finite dimensional tensor for scalar field theories given in ref. [48]. The TRG method, whose computational cost depends on the space-time volume logarithmically, allows us to access to large volume lattices so that the thermodynamic limit and the power law behavior of the susceptibility can be precisely studied. This work is the first study to determine the critical coupling in the continuum limit employing the TRG method.
The dimensionless critical coupling was obtained as λ µ 2 c cont. with sufficiently small error albeit the simplest form of the TRG is employed. The error is the systematic one coming from the finite D effects. Our result shows a reasonable consistency with the recent results obtained by different approaches. The simplest TRG algorithm suffers from the growth of the systematic errors around the criticality. Alternative coarse-graining procedures such as the tensor network renormalization (TNR) [29] and loop-TNR [30] might be useful to obtain more precise results. These methods effectively work around the critical point and, in principle, are applicable to any two dimensional model irrespective of the field contents. The fact that the discretization of scalar fields does not ruin the effectiveness of the TRG approach would encourage further studies using the improved algorithms, and we could expect further improvements for the accuracy of the critical coupling.

A Coarse-graining of tensor network including impurity tensors
We describe the coarse-graining algorithm for a tensor network with an impurity tensor such as Z 1 (K) in eq. (3.13), which is given with a fixed integer D for truncating the SVD.
Before discussing the nonuniform case, we first explain the coarse-graining of a uniform network such as Z (K) in eq. (3.8) without impurity tensors. The graphical representation of a coarse-graining step is given in figure 9. Firstly, by using the SVD, the tensor T with the bond dimension N is decomposed in two ways: (ij)m σ [13] m S [3] m(kl) , (A.1) (li)m σ [24] m S [4] m(jk) , where M [13] and M [24] (N 2 × N 2 matrices) are obtained by arranging the indices of T , σ [13] and σ [24] are the singular values in the descending order, and S [i] (i = 1, 2, 3, 4) denotes the singular matrix. We actually apply eq. (A.1) to the tensors on even sites and eq. (A.2) to ones on odd sites. Figure 9. Coarse-graining step for uniform TN. Circles represent tensors, and closed indices are contracted. A local block on the network is taken into account here. In the first step, every rankfour tensor is decomposed into two rank-three tensors. By contracting the four rank-three tensors, a new rank-four tensor is obtained in the second step. Here, decomposed tensors with tilde are defined byS [1(3)
In the first few steps in the TRG, N 2 might not be larger than D. In this case D is replaced by N 2 for N 2 < D. 12 Note that the coarse-grained tensor T new is defined by gathering S [i] (i = 1, 2, 3, 4) and its bond dimension is at most D. One can see that the tensor network of T is approximately expressed as a coarse network of T new . The number of tensors decreases at each coarse-graining step since T provides two S's while T new is made of four S's. Repeating the procedure above over and over again, we finally find that the network is expressed as a single tensor whose indices are contracted with itself. The value of the tensor network is thus obtained. Now let us turn to the coarse-graining of the tensor network with an impurity tensor. As explained below, an important point is that the number of impurity tensors does not increase beyond four (the impurity tensors are located in four corners of a plaquette) even if the coarse-graining step is repeated many times.
Suppose we consider a coarse-graining procedure for the tensor network with four impurity tensors as shown in figure 10 (left), where the impurity tensors are represented JHEP05(2019)184 Figure 10. Coarse-graining of TN with four impurity tensors. Red, orange, yellow, and green circles denote the impurity tensors, and the normal tensors are shown by blue ones. Note that the impurity tensors colored in red, orange, yellow, and green are inherited in both decomposition and contraction processes. A key point is that the number of impurity tensors does not change through the coarse-graining.
as the red, orange, yellow, and green circles while the normal tensors are given by the blue circles. In the first step, the rank-four tensors are decomposed into rank-three tensors. As seen in the middle panel in figure 10, those rank-three impurity tensors form a closed loop. Thanks to this structure, after taking the contraction, one obtains a network which again contains only four impurity tensors as shown in the right panel in figure 10. In this way, one can suppress the spread of the impurity tensors.
Starting from a tensor network with a single impurity tensor such as eq. (3.13), one obtains a network with two neighboring impurity tensors after a course-graining step. In the next step it is coarse-grained to yield a network with three impurity tensors located on three corners of a plaquette. In this way, the number of impurity tensors stepwisely increases until it reaches four. After that, as shown above, the number of impurity tensors does not change anymore.

B Systematic error from discretization of scalar fields
As shown in section 3 we discretize the scalar fields using the Gauss-Hermite quadrature in order to create a finite dimensional tensor, which gives the approximated expressions for the partition function and expectation values. The order of approximation is characterized by the degree of the Hermite polynomial K. In this appendix we discuss the systematic error associated with the approximation, namely the finite K-effect, for Z and φ . 13 For an one dimensional integral in eq. (3.5), the Gauss-Hermite quadrature gives an exact result when g (x) is a 2K − 1 degree (or less) polynomial function. In the case of multi dimensional integrals, unfortunately, the convergence is not trivial. We have to check whether or not the numerical results do not have large systematic errors from the finite K-effect with K = 256, which is employed in this study.
In figure 11 the K-dependence of φ is presented. It is clear that the results do not depend on K for K ∈ [64, 256] while much larger D-dependence is observed. conclude that our choice K = 256 is large enough, and the main source of systematic errors is the truncation of the SVD for constructing the initial tensor and the coarse-graining steps in the TRG.
Here we also show the D-dependence of ln Z and φ in figures 12 and 13. Although ln Z shows a good convergence, φ is relatively not stable. However, in the large D region (D ∈ [32, 64]), the results seem to be in a fluctuation region, and we set the range of D to [32,64]