Metal-insulator transition in (2+1)-dimensional Hubbard model with tensor renormalization group

We investigate the doping-driven metal-insulator transition of the (2+1)-dimensional Hubbard model in the path-integral formalism with the tensor renormalization group method. We calculate the electron density $\langle n\rangle$ as a function of the chemical potential $\mu$ choosing three values of the Coulomb potential with $U=80$, $8$, and $2$ as representative cases of the strong, intermediate, and weak couplings. We have determined the critical chemical potential at each $U$, where the Hubbard model undergoes the metal-insulator transition from the half-filling plateau with $\langle n\rangle=1$ to the metallic state with $\langle n\rangle>1$. Our results indicate that the model exhibits the metal-insulator transition over the vast region of the finite coupling $U$.

The Hubbard model, which is a simple theoretical model to describe electron systems with repulsive Coulomb interactions, is expected to have rich phase structures so that it has been attracting the interest of not only condensed matter physicists but also particle physicists. It has been widely known that the Hubbard model has a similar path-integral form to the Nambu-Jona-Lasinio (NJL) model [1,2], which is a low energy effective theory in Quantum Chromodynamics (QCD): Both consist of a hopping term and a four-fermi interaction term. Their similarity, unfortunately, leads to sharing the so-called sign problem, which is a notorious difficulty in the numerical analyses based on the Monte Carlo approach.
Recently the authors have successfully applied the tensor renormalization group (TRG) method 1 to investigate the phase transition of the four-dimensional (4d) NJL model at high density and very low temperature [9]. This work was followed by the application of the TRG method to analyze the metal-insulator transition of the (1 + 1)d Hubbard model by calculating the electron density n as a function of the chemical potential µ [12]. Our results for the critical chemical potential µ c and the critical exponent ν are consistent with an exact solution based on the Bethe ansatz [13,14].
In this paper, we apply the TRG method to investigate the doping-driven metal-insulator transition in the (2 + 1)d Hubbard model. 2 The TRG method, which was originally proposed to study two-dimensional (2d) classical spin systems [3], has been developed to study wide varieties of fermionic models in particle physics [5,6,9,[16][17][18][19][20][21][22][23][24]. It is also confirmed that the TRG method does not suffer from the sign problem by studying various quantum field theories [5, 9, 12, 16-18, 20-23, 25-27]. We calculate the electron density n as a function of µ with three choices of U = 80, 8 and 2. The µ dependence of n allows us to determine the critical chemical potential µ c at the doping-driven metal-insulator transition from the halffilling plateau with n = 1 to the metallic state with n > 1. Our results at U = 80, 8 and 2 show that |µ c − U/2| monotonically diminishes as U decreases and seems to converge on |µ c − U/2| = 0 at U = 0. This indicates the possibility that the model exhibits the metal-insulator transition over the wide region of the finite coupling.
This paper is organized as follows. In Sec. 2 we define the Hubbard model in the pathintegral formalism and give a brief description of the numerical algorithm. In Sec. 3 we present the µ dependence of the electron density and determine the critical chemical potential µ c at the doping-driven metal-insulator transition. Section 4 is devoted to summary and outlook.

Formulation and numerical algorithm 2.1. (2+1)-dimensional Hubbard model in the path-integral formalism
We consider the partition function of the Hubbard model in the path-integral formalism on an anisotropic lattice with the physical volume V = L x × L y × β, whose spatial extension is defined as L σ = aN σ (σ = x, y) with a the spatial lattice spacing. β denotes the inverse temperature, which is divided as β = 1/T = ǫN τ . Following Ref. [12], the path-integral expression of the partition function is given by where n = (n x , n y , n τ ) ∈ Λ 2+1 (⊂ Z 3 ) specifies a position on the lattice |Λ 2+1 | = N x × N y × N τ . Since the Hubbard model describes the spin-1/2 fermions, they are labeled by s =↑, ↓, 1 In this paper the TRG method or the TRG approach refers to not only the original numerical algorithm proposed by Levin and Nave [3] but also its extensions [4][5][6][7][8][9][10][11]. 2 The model has also been investigated by the tensor network method based on the Hamiltonian formalism, like a fermionic PEPS, which is also free from the sign problem. For a recent study, see Ref. [15], for example.

2/11
corresponding to the spin-up and spin-down, respectively. Introducing the notation, the action S is defined as The kinetic terms in the spatial directions contain the hopping parameter t. The fourfermi interaction term represents the Coulomb repulsion of electrons at the same lattice site. The chemical potential is denoted by the parameter µ. Note that the half-filling is realized at µ = U/2 in the current definition. We assume the periodic boundary condition in the spatial direction, ψ(N x + 1, n y , n τ ) = ψ(1, n y , n τ ) and ψ(n x , N y + 1, n τ ) = ψ(n x , 1, n τ ), while the anti-periodic one in the temporal direction, ψ(n x , n y , N τ + 1) = −ψ(n x , n y , 1). In the following discussion, we always set a = 1.

Numerical algorithm
Based on Ref. [28], the tensor network representation of Eq. (1) is immediately obtained. Set d = 2 in the Appendix of Ref. [12] and one can find out the Grassmann tensor which generates the Grassmann tensor network of Eq. (1). The resulting Grassmann tensor T ΨxΨyΨτΨτΨyΨx is of rank 6 and we evaluate gTr[ n T ] with the anisotropic TRG (ATRG) algorithm [7] whose extension to the Grassmann integrals, referred as the Grassmann ATRG (GATRG), is given in Ref. [9]. We also follow the coarse-graining procedure employed in the previous study of (1 + 1)d Hubbard model [12]. Firstly, we carry out m τ times of renormalization along with the temporal direction, which can be seen as the imaginary time evolution of the local Grassmann tensor. Secondly, 3d ATRG procedure is applied as the spacetime coarsegraining. As in the case of (1 + 1)d Hubbard model, we have found that the optimal m τ satisfies the condition ǫ2 mτ ∼ O(10 −1 ) in the sense of preserved tensor norm.

Algorithmic-parameter dependence
The partition function of Eq. (1) is evaluated using the numerical algorithm explained above on lattices with the physical volume . We employ U = 80, 8 and 2 for the four-fermi coupling with t = 1 for the hopping parameter. In Fig. 1 on V = 4096 2 × 1677.7216 lattice with ǫ = 10 −4 . In Fig. 2, we plot the D dependence of δ at (U, t) = (8, 1) with the choices of µ = 6.0, 7.5 and 8.5. As we will see below, µ = 6.0 corresponds to n ≈ 1.0 and µ = 8.5 does to n ≈ 1.5. We observe that δ's at these values of µ decrease as a function of D, though some of them are fluctuating.

Strong coupling limit
We first consider the atomic limit at (U, t) = (8, 0). This case is analytically solvable. The electron density n is obtained by the numerical derivative of the thermodynamic potential in terms of µ: In Fig. 3 we compare the numerical and exact results for the µ dependence of n at (U, t) = (8, 0). Note that 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 is at most O(10 −4 ) in the range of 0 ≤ µ ≤ 16. For µ < 0, the exact thermodynamic potential is equal to zero and the results obtained by the TRG are also equal to zero within a level of double precision.
Since the phase diagram of the metal-insulator transition in the (2+1)d Hubbard model is not well known so far, we investigate the µ dependence of n choosing U = 80 as a representative case in the strong coupling region. In Fig. 4 we plot the electron density n as a function of µ in the vicinity of µ ∼ U with D = 80. We have checked that the convergence behavior of δ at U = 80 is better than that at U = 8. We observe that the electron density starts to increase from n = 1 at µ = 77.0(2) or µ/U = 0.9625 (25) and reaches n = 2 with µ 83.0 or µ/U 1.04. The µ dependence of n is smooth and continuous so that there is no signal of the first-order phase transition. We expect that the critical chemical potential µ c at the doping-driven metal-insulator transition approaches to µ c /U = 1 toward the atomic limit and the transition from n = 1 to n = 2 becomes a step-function as a function of µ/U .

Critical chemical potential at U = 8 and 2
Now let us investigate the metal-insulator transition in the intermediate coupling region at (U, t) = (8, 1). There are a lot of previous work to investigate a possible superconducting phase expected in this coupling region. Since we are interested in the thermodynamic and zero-temperature limit, we first check the volume dependence of the electron density n . In Fig. 5 we plot the µ dependence of n at U = 8 changing the lattice sizes with ǫ = 10 −4 , m τ = 12 and D = 80. We observe that the size of (N x , N y , N τ ) = (2 12 , 2 12 , 2 24 ), which corresponds to V = 4096 2 × 1677.7216, is sufficiently large to be identified as the thermodynamic and zero-temperature limit. We observe the n = 0 plateau for µ −4 and the n = 2 one for 12 µ. The half-filling state is characterized by the plateau of n = 1 in the range of 2 µ 6. These plateaus yield the vanishing compressibility κ = ∂ n /∂µ indicating the insulating states. In order to determine the critical chemical potential µ c in the limit of D → ∞ at U = 8 on V = 4096 2 × 1677.7216 lattice we make a global fit of n with D = 80, 72, 64 and 56 in the metallic phase near the transition point. In Fig. 6  degenerate indicating the small D dependence. For the global fit we employ the following quadratic fitting function: with µ c (D) = µ c (D = ∞) + γ/D, where α, β, γ and µ c (D = ∞) are the fit parameters. The solid curves in Fig. 6 represent the fit results over the range of 6.3 ≤ µ ≤ 8.0. We obtain µ c (D = ∞) = 6.43(4), which is presented in Table 1 together with other fitting results. It may be instructive to compare Figs. 5, 6, and the estimated location of µ c (D = ∞) with numerical data in Refs. [29,30], though their calculations are carried out on very small lattice sizes and at low but finite temperatures. We repeat the same analysis for the weak coupling case at U = 2. We apply the fit function of Eq. (7) to four data sets with the bond dimensions of D = 80, 72, 64 and 56. Fit results are depicted in Fig. 7 and their numerical values are summarized in Table 1. Our results show that the deviation of |µ c (D) − U/2| diminishes as the Coulomb potential U decreases. It is likely that |µ c (D) − U/2| vanishes only at U = 0. This means that the model exhibits the metal-insulator transition over the wide regime of the finite coupling, including the weak coupling region. This conclusion may provide us a different scenario of the phase diagram from that predicted by the dynamical mean-field theory (DMFT) [31]; there exists some U c such that no metal-insulator transition occurs with U < U c .

Summary and outlook
We have investigated the doping-driven metal-insulator transition of the (2+1)d Hubbard model in the path-integral formalism employing the TRG method. The electron density n is calculated in the wide range of µ corresponding to 0 ≤ n ≤ 2. We have also determined the critical chemical potential µ c at three values of U . Our results indicate that the deviation |µ c − U/2| vanishes only at U = 0. This means that the model exhibits the metal-insulator transition over the vast regime of the finite coupling U . As a next step, it would be interesting to investigate the metal-insulator transition of the (3+1)d Hubbard model.