Tensor renormalization group approach to (1+1)-dimensional Hubbard model

We investigate the metal-insulator transition of the (1+1)-dimensional Hubbard model in the path-integral formalism with the tensor renormalization group method. The critical chemical potential $\mu_{\rm c}$ and the critical exponent $\nu$ are determined from the $\mu$ dependence of the electron density in the thermodynamic limit. Our results for $\mu_{\rm c}$ and $\nu$ show consistency with an exact solution based on the Bethe ansatz. Our encouraging results indicate the applicability of the tensor renormalization group method to the analysis of higher-dimensional Hubbard models.


I. INTRODUCTION
The tensor renormalization group (TRG) method 1 , which was originally proposed to study two-dimensional (2d) classical spin systems in the field of condensed matter physics [1], has been now used to study wide varieties of models in particle physics taking several advantages over the Monte Carlo method. (i) The TRG method does not suffer from the sign problem as already confirmed by studying various quantum field theories [3,[7][8][9][10][11][12][13][14]. (ii) Its computational cost depends on the system size only logarithmically. (iii) It allows direct manipulation of the Grassmann variables [3,4,7,15]. (iv) We can obtain the partition function or the path-integral itself.
The sign problem is common both in particle physics and condensed matter physics. A typical example in particle physics is the lattice QCD at finite density, where an introduction of the chemical potential causes the sign problem, and the Hubbard model is notorious in the condensed matter physics. Recently the authors have successfully applied the TRG method to analyze the phase transition of the 4d Nambu−Jona-Lasinio (NJL) model at high density and very low temperature [7]. The study of the NJL model has two important aspects. Firstly, the NJL model is a prototype of QCD. Their phase structures are expected to be similar so that the study of the NJL model at finite density is a good testbed before investigating the finite density QCD. Secondly, the NJL model has a similar path-integral form to the Hubbard model: Both consist of a hopping term and a four-fermi interaction term. This indicates that the technical details of the TRG method employed in the analysis of the NJL model could apply to the Hubbard model. It is interesting to investigate whether or not the TRG method overcomes the sign problem in the Hubbard model.
In this paper, we investigate the metal-insulator transition of the (1+1)d Hubbard model in the path-integral formalism by calculating the electron density as a function of the chemical potential µ. After examining the imaginary-time discretization effects and the temperature dependence, we determine the critical value of the chemical potential µ c and the critical exponent ν in the thermodynamic limit at the zero temperature. Our results for µ c and ν show agreement with the theoretical prediction based on the Bethe ansatz [16,17].
This paper is organized as follows. In Sec. II we define the Hubbard model in the path-integral formalism and explain the numerical algorithm. In Sec. III we show our results and compare them with theoretical predictions. Section IV is devoted to summary and outlook.
Since the model describes spin-1/2 particles, the spatial auxiliary Grassmann field Ψ σ has 2 (hopping terms) × 2 (spin d.o.f.) components and the temporal one Ψ τ has 1 × 2 components. Using the Grassmann tensor T in Eq (4), the path integral Z is expressed by See Appendix A for the detailed explanation to derive the above Grassmann tensor and its tensor network.

C. Numerical algorithm
We employ the higher-order TRG (HOTRG) algorithm [2] to evaluate the Grassmann tensor network in Eq (5). Using the HOTRG, we firstly carry out m τ times of renormalization along the temporal direction. This procedure converts the initial Grassmann tensor T ΨσΨτΨτΨσ into the coarse-grained one T ΞσΨτΨτΞσ . Secondly, we em-ploy the 2d HOTRG procedure, regarding T ΞσΨτΨτΞσ as the initial tensor, to obtain the coarse-grained Grassmann tensor T Ξ σ Ψ τΨ τΞ σ . Note that with sufficiently small (< 1), little truncation error is accumulated with the first m τ times of renormalization along τ -direction. This is because the contribution from the spatial hopping terms, which are of O( ), is smaller than that from the temporal one, which is of O(1). For the (1 + 1)d Hubbard model, we found that the optimal m τ satisfied the condition 2 mτ ∼ O(10 −1 ). 3 When one applies the TRG approach to evaluate the path integral over the Grassmann fields, it is practically useful to encode the Grassmann parity of the auxiliary Grassmann fields into the subscripts of the coefficient tensor. We identify the coefficient tensor in Eq. (4) as a four-rank tensor T xtx t , where x, x = 1, · · · , 2 4 and t, t = 1, · · · , 2 2 . These new indices are defined as in Tables II C and II C. Notice that x(x ) = 1, · · · , 8 correspond to the Grassmann-even sector and x(x ) = 9, · · · , 16 the Grassmann-odd one in Ψ σ (Ψ σ ). Similarly, t(t ) = 1, 2 correspond to the Grassmann-even sector and t(t ) = 3, 4 the Grassmann-odd one in Ψ τ (Ψ τ ). These mappings help us to carry out the singular value decompositions with some block-diagonal representations as explained in Ref. [7].

III. NUMERICAL RESULTS
The partition function of Eq. (1) is evaluated using the numerical algorithm explained above on lattices with the with the periodic boundary condition for the spacial direction and the anti-periodic one for the temporal direction. We employ t = 1 for the hopping parameter and U = 4 for the four-fermi coupling. In Fig. 1 we plot the µ dependence of the thermodynamic potential ln Z/V on V = L × β = 4096 × 1677.7216 with the bond dimension D = 80 in the HOTRG algorithm choosing = 2 12 × 10 −4 , 2 8 × 10 −4 , 2 4 × 10 −4 , 10 −4 . For each value of , m τ is decided via the condition 2 mτ = 2 12 ×10 −4 = O(10 −1 ). We find clear discretization effects 3 A similar remark is also mentioned in Ref. [2], where the 3d HOTRG is applied to 2d quantum transverse Ising model in the path-integral formalism.
for the = 2 12 ×10 −4 case. On the other hand, the results with = 2 4 ×10 −4 and 10 −4 show good consistency. This means that the discretization effects with = 10 −4 are negligible. We investigate the convergence behavior of the thermodynamic potential defining the quantity on V = 4096 × 1677.7216 lattice with = 10 −4 . In Fig. 2, we plot the D dependence of δ at µ = 2.75 and 2.00, which are near and far away from the critical point µ c , respectively, as we will see below. We observe that  Before presenting the U/t = 4 results let us consider the (U, t) = (4, 0) and (0, 1) cases. Since these cases are analytically solvable, it is instructive to compare the numerical results for the electron density with the exact ones. The electron density n is obtained by the numerical derivative of the thermodynamic potential in terms of µ: In Figs. 3 and 4 we compare the numerical and analytic results for the µ dependence of n . In both cases we observe good consistencies over the wide range of µ.
Note that for the case of (U, t) = (4, 0) in Fig. 3, we set m τ = 24 because this case is equivalent to the model defined on V = 1 × β lattice. Thanks to the vanishing hopping structure in the spatial direction, we can always perform an exact tensor contraction in Eq. (5). In Fig. 4 we employ finer resolution of µ around 1 |µ| 2 in order to follow the complicated µ dependence of n .  The results indicate that the size (N σ , N τ ) = (2 12 , 2 24 ), which corresponds to V = 4096 × 1677.7216, is sufficiently large to be identified as the thermodynamic and zero-temperature limit.
The half-filling state is characterized by the plateau with n = 1 in the range of 1.3 µ 2.7. We also observe the continuous change from n = 1 to n = 2 over the range of 2.7 µ 6.5. Figure 6 shows µ dependence of n near the criticality on V = 4096 × 1677.7216. The abrupt change of n at µ ≈ 2.70 in Fig. 6 indicates a metal-insulator transition.   We determine the critical chemical potential µ c (D) and the critical exponent ν on V = 4096 × 1677.7216 lattice by fitting n in the metallic phase around the transition point with the following form: where A, B, µ c (D) and ν are the fit parameters. The solid curve in Fig. 6 shows the fitting result over the range of 2.68 ≤ µ ≤ 3.00. We obtain µ c (D) = 2.698(1) and ν = 0.51(2) at D = 80. Our result for the critical exponent is consistent with the theoretical prediction of ν = 1/2. A previous Quantum Monte Carlo simulation with small spatial extension up to L = 24 also yielded the same conclusion [22]. In order to extrapolate the result of µ c (D) to the limit D → ∞, we repeat the calculation changing D. The  Fig. 7

IV. SUMMARY AND OUTLOOK
We have investigated the metal-insulator transition of the (1+1)d Hubbard model in the path-integral formalism employing the TRG method. Extrapolating µ c (D) to the limit D → ∞, we have estimated the critical chemical potential, which shows good consistency with the theoretical prediction based on the Bethe ansatz. We have determined the critical exponent ν, which is also consistent with the exact solution. These encouraging results show the effectiveness of the TRG approach for the study of the Hubbard model and the related fermion models being free from the sign problem. It is worth emphasizing that the TRG approach is efficient not only in the lowerdimensional systems but also in the higher-dimensional ones, as confirmed in the earlier works [2, 4-7, 14, 15, 23-26]. As a next step, we are planning to investigate the phase diagram of the higher-dimensional Hubbard models, improving the TRG method successfully applied in this work.