On the Kinetic Energy Driven Superconductivity in the Two-Dimensional Hubbard Model

We investigate the role of kinetic energy for the stability of superconducting state in the two-dimensional Hubbard model on the basis of an optimization variational Monte Carlo method. The wave function is optimized by multiplying by correlation operators of site off-diagonal type. This wave function is written in an exponential-type form given as ψλ = exp(−λK)ψG for the Gutzwiller wave function ψG and a kinetic operator K. The kinetic correlation operator exp(−λK) plays an important role in the emergence of superconductivity in large-U region of the two-dimensional Hubbard model, where U is the on-site Coulomb repulsive interaction. We show that the superconducting condensation energy mainly originates from the kinetic energy in the strongly correlated region. This may indicate a possibility of high-temperature superconductivity due to the kinetic energy effect in correlated electron systems.


Introduction
It is important and challenging to clarify the mechanism of high-temperature superconductivity in cuprates. It has been studied intensively for more than three decades [1]. In the study of cuprate high-temperature superconductivity, it is important to understand the ground state phase diagram. For this purpose, we should investigate electronic models with strong correlation.
We employ an optimization variational Monte Carlo method to investigate the ground state of the 2D Hubbard model. In a variational Monte Carlo method, we use variational wave functions with strong correlation between electrons [33,34,[36][37][38][39]. A variational wave function is improved by introducing correlation operators. The wave function used in this paper is obtained by multiplying the Gutzwiller function by exp(−λK) operators, where K is the kinetic part of the Hamiltonian [47,48,63,64]. We can optimize the wave function further by multiplying by exponential-type operators again [63]. The ground-state energy is lowered greatly with this wave function [47].
The kinetic energy-driven mechanism of superconductivity has been examined in the study of cuprate superconductors for the Hubbard model [42,[65][66][67] and the t-J model [68][69][70]. This is an interesting subject since there is a possibility that kinetic energy pairing occurs in high-temperature cuprates. In this paper, we discuss the role of kinetic energy based on the improved wave function for the Hubbard model. We evaluate the superconducting condensation energy as a sum of the kinetic and Coulomb energy contributions, and we show that the kinetic energy contribution dominates the condensation energy in the strongly correlated region of large U.

Hubbard Hamiltonian
We investigate the two-dimensional Hubbard model. The Hamiltonian of the Hubbard model is where t ij indicates the transfer integral, and U is the strength of the on-site Coulomb interaction. The transfer integral is t ij = −t when i and j are nearest-neighbor pairs ij . We put t = 0 in this paper, where t ij = −t when i and j are next-nearest neighbor pairs. N and N e denote the number of lattice sites and the number of electrons, respectively. The energy unit is given by t. We define the non-interacting part as K: When two electrons with spin up and down are at the same site, the energy becomes higher by U. This simple interaction may cause many interesting phenomena, such as the metal-insulator transition, antiferromagnetic magnetism, and superconductivity. The metal-insulator transition occurs at half filling when U is as large as the bandwidth. The effective Hamiltonian is given by the Heisenberg model when U is large in the halffilled case, which leads to the t-J model when holes are doped. Magnetic properties of materials may be described by the Hubbard model by introducing suitable magnetic orders. We discuss superconductivity in the strongly correlated region in this paper. The emergence of superconducting state in this region is closely related to the kinetic energy gain that increases as U increases.

Optimized Wave Function
In a variational Monte Carlo method, we calculate the expectation values of physical properties by using a Monte Carlo procedure. We start with the Gutzwiller wave function to take account of electron correlation. The Gutzwiller wave function is where P G is the Gutzwiller operator P G = ∏ j (1 − (1 − g)n j↑ n j↓ ), where g is the variational parameter in the range of 0 ≤ g ≤ 1. ψ 0 indicates a one-particle state, such as the Fermi sea, the BCS state, and the antiferromagnetic state. We improve the wave function by multiplying by the correlation operator given as [47,63,[71][72][73][74][75] where K is the kinetic part of the Hamiltonian, and λ is a real variational operator [39,63,72].
There are other methods to improve the Gutzwiller wave function by using Jastrowtype operators [41,76,77]. This operator is written as where d j is the operator for the doubly-occupied site given as d j = n j↑ n j↓ , and e j is that for the empty site given by e j = (1 − n j↑ )(1 − n j↓ ). η is the variational parameter that takes the value in the range of 0 ≤ η ≤ 1. τ indicates a vector connecting nearest-neighbor sites. The Jastrow-type wave function is written as ψ J = P Jdh ψ G . An important difference between ψ λ and ψ J is that P Jdh is the site-diagonal operator, while ψ λ is the site-off diagonal operator. The one-particle state ψ 0 is written in the form where {ϕ 0 j } is a set of basis functions of the one-particle state in the site representation on a lattice, with j being the label for the electron configuration. ψ λ and ψ J are given as Since P G and P Jdh are diagonal operators, ψ J is written as where a J j = a J j (g, η) is determined by parameters g and η. ψ J is given by a wave function, where the coefficients {a 0 j } are modified in the one-particle state. Instead, ψ λ is not so simple because e −λK generates other basis states from ϕ 0 j , which means that off-diagonal elements ϕ 0 i Kϕ 0 j are effectively taken into account. We write ψ λ as where the set of basis states {ϕ j } may contain basis states which are not included in {ϕ 0 j } since some coefficients a 0 s may vanish accidentally.
The expectation values are calculated numerically by using the auxiliary field method following a Monte Carlo algorithm [63,78]. We show the ground-state energy per site as a function of U for ψ G and ψ λ in Figure 1. The energy is lowered due to the exponential factor e −λK . In the region of large U, the energy lowering mainly comes from the kinetic energy gain, as shown later.

Correlated Superconducting State
The correlated superconducting state is obtained by multiplying the BCS wave function by a correlation operator. The BCS wave function is with coefficients u k and v k that appear in the ratio u k /v k = ∆ k /(ξ k + ξ 2 k + ∆ 2 k ), where ∆ k is the gap function with k dependence, and ξ k = k − µ is the dispersion relation of conduction electrons. We adopt the d-wave symmetry for ∆ k , and we do not consider other symmetries in this paper since an SC state with other symmetry unlikely becomes stable in the simple single-band Hubbard model. The Gutzwiller BCS state is given by where P N e indicates the operator to extract the state with N e electrons. Here, the total electron number is fixed, and the chemical potential in ξ k is regarded as a variational parameter.
The improved superconducting wave function is written as In the formulation of ψ λ , we cannot fix the total electron number, and we should use a different Monte Carlo sampling procedure. We perform the electron-hole transformation for down-spin electrons: and not for up-spin electrons: c k = c k↑ . In the real space, we have c i = c i↑ and d i = c † i↓ . The electron pair operator c † k↑ c † −k↓ is transformed to the hybridization operator c † k d k . Then, we can use the auxiliary field method after this transformation in a Monte Carlo simulation. This is performed by expressing the Gutzwiller operator in the form [72].
where g = e −α , cosh(2a) = e −α/2 , and s i is the auxiliary field that takes the value of ±1. By using this form, the expectation value is calculated as a sum of terms with respect to auxiliary fields, for which we apply the Monte Carlo procedure [63,72].

Kinetic Energy in the Superconducting State
In this section, we discuss the role of kinetic energy in the superconducting state in the two-dimensional Hubbard model. We consider the large-U region, where the kinetic energy of electrons would play an important role. Here, large-U means that U is much larger than the band width. The ground-state energy is determined by the balance of the kinetic energy and the Coulomb energy. We show the energy expectation values as a function of λ in Figure 2, where E g /N is the ground-state energy per site, E kin /N is the expectation value of the non-interacting part of the Hamiltonian E kin = ψ λ Kψ λ / ψ λ ψ λ , and E U denotes the Coulomb energy given by E U = U ψ λ ∑ i n i↑ n i↓ ψ λ / ψ λ ψ λ . E g has a minimum value for a finite value of λ.
The kinetic energy part gives a large contribution to E g when U is large. This is shown in Figure 3. When U > 10t, E U for ψ λ almost agrees with that for ψ G . The difference of E kin for ψ λ and ψ G increases for U > 10t. The region that may be called the 'kinetic energy phase' exists when 10 < U/t. Let us consider the the difference of the kinetic energy defined as where E kin (ψ G ) and E kin (ψ λ ) denote the kinetic energy for ψ G and ψ λ , respectively. Since ψ G = ψ λ=0 with vanishing λ, we can write ∆E kin = E kin (λ = 0) − E kin (λ) for the optimized value λ. We show ∆E kin /N as a function of U in Figure 4, where the hole doping rate x is x = 0.12. ∆E kin begins to increase when U is of the order of the band width U ∼ 8t. ∆E kin has a broad peak when 15t < U < 20t. We define the SC condensation energy ∆E sc and the kinetic energy condensation energy ∆E kin−sc as where ∆ = ∆ sc is the superconducting order parameter, and ∆ opt is its optimized value. We assumed the d-wave symmetry for ∆ k : ∆ k = ∆(cos k x − cos k y ). ∆E kin−sc /N is also shown in Figure 4. The figure indicates that ∆E kin−sc increases for large U, showing a similar behavior to ∆E kin . ∆E kin−sc can change the sign when U is small, which is consistent with the analysis for Bi 2 Sr 2 CaCu 2 O 8+δ [79].
In Figure 5, we show ∆E kin /N and ∆E kin−sc /N for x = 0.20. ∆E kin for x = 0.20 is smaller than that for x = 0.12. The kinetic energy condensation energy ∆E kin−sc is also reduced for x = 0.20 compared to that for x = 0.12. This result inevitably leads to the decrease of the SC condensation energy [48]. Hence, the kinetic energy effect becomes weak when the hole doping rate is large.
The Coulomb energy gain in the presence of the SC order parameter ∆E U−sc is also evaluated, where ∆E U−sc is the Coulomb energy contribution to the SC condensation energy; namely, we have Our result shows that for the improved state ψ λ−BCS . We show the SC condensation energy ∆E sc (∆ sc ) and the Coulomb energy part ∆E U−sc (∆ sc ) as a function of ∆ sc in Figure 6 where ∆E sc (∆) = E g (∆ = 0) − E g (∆) and ∆E U−sc (∆) = E U (∆ = 0) − E U (∆). The positive ∆E kin−sc > 0 indicates that the SC state ψ λ−BCS becomes stable due to the kinetic energy effect.

Summary
We investigated the ground state of the two-dimensional Hubbard model by using the optimization variational Monte Carlo method with focus on the strongly correlated large-U region. The ground-state energy is greatly lowered due to the exp(−λK) correlation operator. The optimization variational Monte Carlo method is effective even in strongly correlated regions where U is much larger than the bandwidth.
In the large-U region, the kinetic energy term becomes dominant in the ground state, and the kinetic energy effect dominates the antiferromagnetic correlation. The kinetic energy difference ∆E kin = E kin (ψ G ) − E kin (ψ λ ) increases when U > 10t and has a broad peak. The kinetic condensation energy ∆E kin−sc behaves like the kinetic energy difference ∆E kin for U > 10t. There is a correlation between ∆E kin−sc and ∆E kin . The condensation energy ∆E sc mainly comes from the kinetic energy part ∆E kin−sc and there is a competition between the kinetic energy and the Coulomb energy. Hence there are two competitions in the Hubbard model; one is that between superconductivity and antiferromagnetism and the other is that between kinetic energy effect and Coulomb repulsive interaction. As a result of competitions, superconducting transition occurs. The result shows that superconductivity in the strongly correlated region is induced by kinetic energy effect. We expect that high-temperature superconductivity would be realized in the strongly correlated region of the two-dimensional Hubbard model. Acknowledgments: A part of computations was supported by the Supercomputer Center of the Institute for Solid State Physics, the University of Tokyo.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: VMC variational Monte Carlo method AF antiferromagnetic SC superconductivity or superconducting 2D two-dimensional AFI antiferromagnetic insulator PI paramagnetic insulator