Antiferromagnetism, Superconductivity and Phase Diagram in the Two-Dimensional Hubbard Model -- Off-Diagonal Wave Function Monte Carlo Studies of Hubbard Model III --

TL;DR

This study uses an advanced variational Monte Carlo method to explore the phase diagram of the 2D Hubbard model, revealing the conditions favoring antiferromagnetism and superconductivity, with off-diagonal wave functions playing a key role.

Contribution

It introduces an off-diagonal wave function approach in variational Monte Carlo to analyze competing phases in the Hubbard model, highlighting the importance of many-body effects.

Findings

Superconductivity is most favored at zero next-nearest neighbor hopping ($t'=0$).

The phase diagram includes both antiferromagnetic and superconducting phases.

Many-body effects significantly influence spin correlations and superconductivity.

Abstract

We investigate the ground-state phase diagram of the two-dimensional Hubbard model based on the optimization variational Monte Carlo method. We use a wave function that is an off-diagonal type given as $\psi=\exp(-\lambda K)P_G\psi_0$, where $\psi_0$ is a one-particle state, $P_G$ is the Gutzwiller operator, $K$ is the kinetic operator, and $\lambda$ is a variational parameter. The many-body effect plays an important role as an origin of spin correlation and superconductivity in correlated electron systems. We examine the competition between the antiferromagnetic state and superconducting state by varying the Coulomb repulsion $U$, the band parameter $t'$ and the electron density $n_e$. We show a phase diagram that includes superconducting and antiferromagnetic phases and that $t'=0$ is most favorable for superconductivity.

Equations

$$H=\sum _{ij\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$${\psi }_G=P_G{\psi }_0,$$

$${\psi }_{\lambda }=\exp (-\lambda K){\psi }_G,$$

$$P_{Jdh}=\prod_{j}\left(1-(1-\eta)\prod_{\tau}\Big{[}d_{j}(1-e_{j+\tau})+e_{j}(1-d_{j+\tau})\Big{]}\right),$$

$${\psi }_{\eta }=P_{Jdh}{\psi }_G.$$

$$H_{AF}=\sum _{ij\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }-{\Delta }_{AF}\sum _{i\sigma }(-1{)}^{x_i+y_i}\sigma n_{i\sigma },$$

$${\psi }_{\lambda ,AF}=\exp (-\lambda K)P_G{\psi }_{AF}.$$

$${\psi }_{BCS}=\prod _k(u_k+v_k{c}_{k\uparrow }^{†}{c}_{-k\downarrow }^{†})|0⟩,$$

$${\psi }_{G-BCS}=P_{N_e}{P}_G{\psi }_{BCS},$$

$${\psi }_{\lambda }=e^{-\lambda K}{P}_G{\psi }_{BCS}.$$

$$d_k=c_{-k\downarrow }^{†},d_k^{†}=c_{-k\downarrow },$$

$$\frac{1}{{\chi }_c}=\frac{{\partial }^{2}E(N_e)}{\partial N_e^2}=\frac{E(N_e+\delta N_e)+E(N_e-\delta N_e)-2E(N_e)}{(\delta N_e{)}^2},$$

$${\delta }^{2}E(N_e)\equiv E(N_e+\delta N_e)-2E(N_e)+E(N_e-\delta N_e)$$

$$\Delta E=E({\Delta }_{SC}=0)-E({\Delta }_{SC,\mathrm{opt}}),$$

$${\psi }_{\lambda }=e^{-\lambda K}{P}_G{\psi }_{BCS-AF}.$$

$${\psi }^{(4)}=e^{-{\lambda }^{\prime }K}{P}_G(g^{\prime })e^{-\lambda K}{P}_G(g){\psi }_{AF},$$

Full text

1 Introduction

The mechanism and properties of high-temperature superconductivity have been studied vigorously for more than 30 years since the discovery of cuprate high-temperature superconductors[1]. High-temperature cuprates are typical strongly correlated systems since the parent materials are Mott insulators when no carriers are doped. It is important to understand the electronic properties of strongly correlated electron systems because high-temperature cuprates are typical strongly correlated systems and the parent materials are Mott insulators when no carriers are doped.

The CuO2 plane commonly contained in high-temperature cuprates consists of oxygen atoms and copper atoms. The electronic model for this plane is given by the d-p model or three-band Hubbard model[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. It appears very difficult to understand the ground-state phase diagram of the d-p model because of strong correlation between electrons. The two-dimensional (2D) single-band Hubbard model[25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45] has been investigated as a simplified model of the d-p model. The ladder model[46, 47, 48, 49, 50, 51] has also been studied in relation to the mechanism of superconductivity in a correlated electron system.

The Hubbard model is one of the fundamental models in condensed matter physics. It was first introduced to understand the metal-insulator transition[52] and is also used to describe the magnetic properties of various compounds[53, 54]. By employing the Hubbard model, we can understand the appearance of inhomogeneous states such as stripes[55, 56, 57, 58, 59, 60, 61, 62] and checkerboard-like density wave states[63, 64, 65, 66], whose existence was reported for high-temperature cuprates.

Recent studies on the 2D Hubbard model indicate that a superconducting (SC) phase exists in the ground state[45]. This shows the possibility that the 2D Hubbard model may account for high-temperature superconductivity. We show the order parameters of the antiferromagnetic (AF) state and SC state as a function of the interaction parameter $U$ in Fig. 1. The result shows that high-temperature superconductivity may occur in the strongly correlated region of the Hubbard model where the interaction $U$ is greater than the bandwidth.

A variational Monte Carlo method is a useful tool to investigate the electronic properties of strongly correlated electron systems when we calculate the expectation values numerically[32, 33, 34, 35, 36, 37]. In general, a variational wave function is improved by introducing new variational parameters to control the electron correlation. In our method the wave functions are optimized by multiplying an initial function by $\exp(-S)$-type operators[45, 67], where $S$ is a suitable correlation operator. The Gutzwiller function is also written in this form. An optimization process is performed in a systematic way by multiplying by the exponential-type operators repeatedly[67]. The ground-state energy is indeed lowered considerably by using this type of wave function[45].

In this paper we investigate the stability of the antiferromagnetically ordered state and show the phase diagram of the ground state of the 2D Hubbard model. In the strongly correlated region, the AF correlation is suppressed and the SC correlation is enhanced. Near the boundary of the AF region, a large spin fluctuation is induced, which is considered to give rise to high-temperature superconductivity. The paper is organized as follows. In Sect. 2, we present the model Hamiltonian and wave functions that we use in the optimization variational Monte Carlo method. In Sect. 3, we examine the wave function with the AF order parameter to show the region where the AF state is stabilized. In Sect. 4, we discuss the phase separation that may occur near half-filling. In Sect. 5, we investigate the SC state and show the phase diagram. We give a summary in the final section.

2 Optimization Variational Monte Carlo Method

2.1 Hamiltonian

The Hubbard model is written as

[TABLE]

where $t_{ij}$ indicates the transfer integral and $U$ is the strength of the on-site Coulomb interaction. We set $t_{ij}=-t$ when $i$ and $j$ are nearest-neighbor pairs $\langle ij\rangle$ and $t_{ij}=-t^{\prime}$ when $i$ and $j$ are next-nearest-neighbor pairs. We consider this model in two dimensions, and $N$ and $N_{e}$ denote the number of lattice sites and the number of electrons, respectively. The energy unit is given by $t$.

2.2 Off-diagonal wave function

In a variational Monte Carlo method, we employ a wave function that is suitable for the system being considered and evaluate the expectation values by using a Monte Carlo procedure. To take into account the correlation between electrons, we start from the Gutzwiller wave function given by

[TABLE]

where $P_{G}$ is the Gutzwiller operator $P_{G}=\prod_{j}(1-(1-g)n_{j\uparrow}n_{j\downarrow})$, where $g$ is the variational parameter in the range of $0\leq g\leq 1$. $\psi_{0}$ indicates a trial one-particle state.

Because the Gutzwiller function is very simple and is not enough to take account of electron correlation, we should improve the wave function. There are several methods to improve the wave function. One method is to multiply the Gutzwiller function by an exponential-type operator. The wave function is written as[45, 67, 68, 69, 70, 71, 72, 73]

[TABLE]

where $K$ is the kinetic part of the Hamiltonian and $\lambda$ is a real variational operator[37, 67, 69]. The expectation values are calculated by using the auxiliary field method[67, 74]. The other method is to introduce a Jastrow-type operator[39]. We control the nearest-neighbor correlation by multiplying by the operator

[TABLE]

where $d_{j}$ is the operator for the doubly occupied site given as $d_{j}=n_{j\uparrow}n_{j\downarrow}$ and $e_{j}$ is that for the empty site given by $e_{j}=(1-n_{j\uparrow})(1-n_{j\downarrow})$. $\eta$ is the variational parameter in the range $0\leq\eta\leq 1$. With this operator, we can include the doublon-holon correlation:

[TABLE]

It is possible to generalize the Jastrow operator to consider long-range electron correlation by introducing new variational parameters[75, 76].

In this paper we use the wave function of the exponential type in eq. (3). We call this type of wave function the off-diagonal wave function since the off-diagonal correlation in the site representation is taken into account in this wave function. We believe that it is more important to consider off-diagonal electron correlation than the diagonal electron correlation. In fact, the energy is further lowered when we employ the off-diagonal wave function[45].

2.3 Antiferromagnetic state

The AF one-particle state $\psi_{AF}$ is given by the eigenfunction of the AF trial Hamiltonian

[TABLE]

where $\Delta_{AF}$ is the AF order parameter and $(x_{i},y_{i})$ represents the coordinates of site $i$. The wave function is written as

[TABLE]

In general, the AF state is very stable in the Hubbard model near half-filling. Thus, it is important to control the AF magnetic order so that the SC state is stabilized and realized.

The stability of the AF state depends mainly on the electron density $n_{e}$, the interaction strength $U$, the transfer integral $t^{\prime}$, and long-range transfers in the single-band Hubbard model. The AF correlation is induced as $U$ increases from zero in the weakly correlated region and is maximized when $U$ is of the order of the bandwidth, say at $U=U_{c}$, when carriers are doped. When $U$ becomes larger than $U_{c}$, the AF correlation starts to decrease. In the region where $U$ is extremely large, the AF correlation is suppressed to a small value by the large fluctuation. This is shown in Fig. 1. This is a crossover between the weakly correlated region and strongly correlated region.

2.4 Correlated superconducting state

The SC state is represented by the BCS wave function

[TABLE]

with coefficients $u_{k}$ and $v_{k}$ that appear in the ratio $u_{k}/v_{k}=\Delta_{k}/(\xi_{k}+\sqrt{\xi_{k}^{2}+\Delta_{k}^{2}})$, where $\Delta_{k}$ is the gap function with ${\bf k}$ dependence and $\xi_{k}=\epsilon_{k}-\mu$ is the dispersion relation of conduction electrons. We assume $d$-wave symmetry for $\Delta_{k}$: $\Delta_{k}=\Delta_{SC}(\cos k_{x}-\cos k_{y})$. The Gutzwiller BCS state is formulated as

[TABLE]

where $P_{N_{e}}$ indicates the operator used to extract the state with $N_{e}$ electrons. In this wave function the electron number is fixed and thus the chemical potential in $\xi_{k}$ is regarded as a variational parameter. In the formulation of $\psi_{\lambda}$, we use the BCS wave function without fixing the total electron number, namely, without the operator $P_{N_{e}}$. The chemical potential $\mu$ in $\xi_{k}$ is not a variational parameter and is used to adjust the total electron number. The wave function is given as

[TABLE]

We perform the electron-hole transformation for down-spin electrons:

[TABLE]

and not for up-spin electrons: $c_{k}=c_{k\uparrow}$. The electron pair operator $c^{{\dagger}}_{k\uparrow}c^{{\dagger}}_{-k\downarrow}$ denotes the hybridization operator $c^{{\dagger}}_{k}d_{k}$ in this formulation.

3 Antiferromagnetic phase

Since the SC state competes with the AF state, it is important to clarify the region of the AF state in the parameter space. We use $U$ and $t^{\prime}$ to control the strength of the AF correlation. The Coulomb interaction $U$ is important since the magnitude of the AF magnetism can be controlled by changing $U$. The transfer integral $t^{\prime}$ is also important and shows nontrivial effect on the stability of the AF magnetic order. One may expect that the AF region will be small when including $t^{\prime}$ in the model. This is not, however, true. As $|t^{\prime}|$ increases, the AF correlation increases, where we assume negative $t^{\prime}$ in this paper. From the viewpoint of competition between superconductivity and AF ordering, $t^{\prime}=0$ is most favorable for superconductivity.

We show the condensation energy $\Delta E_{AF}$ due to the AF magnetic order as a function of $1-n_{e}$ in Fig. 2 for $t^{\prime}=0$. When $U$ is as large as $U/t\geq 14$, the AF region exists up to 10$\%$ doping. When $t^{\prime}=-0.2t$, the AF region expands up to about 20$\%$ doping where $1-n_{e}\sim 0.2$ even for large $U$. We show this in Fig. 3. The AF region becomes larger as $|t^{\prime}|$ increases. Figure 4 shows $\Delta E_{AF}$ as a function of $x$ for $t^{\prime}=0$, $-0.1$, and $-0.2$, where we use $U/t=18$. We show the AF region on the $x-t^{\prime}$ plane in Fig. 5. The AF state dominates near half-filling and is stabilized as $|t^{\prime}|$ increases. The $d$-wave SC state exists near the boundary in Fig. 5.

4 Phase Separation

We discuss the phase separation in the 2D Hubbard model in this section. The existence of phase separation near half-filling in the 2D Hubbard model has been discussed[77, 78, 79]. This is a subject concerning the AF correlation and charge distribution in the case of small doped carriers. When there is a strong AF correlation between neighboring electrons, there may be a tendency that doped holes form clusters due to an effective attractive interaction between electrons. This means the possibility of phase separation with clusters of holes, depending on the AF correlation, attractive interaction, and kinetic energy gain. This is similar to the instability in the t-J model[80].

We examine an instability toward phase separation by evaluating the charge susceptibility

[TABLE]

where $E(N_{e})$ is the ground-state energy and $N_{e}$ is the number of electrons. This is proportional to the second derivative of the energy $E(N_{e})$ with respect to the electron number. The negative sign of $\chi_{c}$ indicates an instability toward the phase separation. This instability is very subtle. Once the phase separation occurs, the ground state becomes an insulating state.

We show the energy as a function of the doping rate $x$ for $t^{\prime}=0$ in Fig. 6. The curve of the energy is usually convex downward, that is, $\chi_{c}>0$, but the sign of $\chi_{c}$ changes in the region near half-filling. We show

[TABLE]

for $N_{e}=2$ in Fig. 7. This indicates that there is an instability toward the phase separation when $x<0.06$ for $U/t=18$ and $t^{\prime}=0$. This is similar for $U/t=14$ and 12. When $t^{\prime}$ is nonzero and negative, the instability toward the phase separation occurs for a smaller doping rate. We show $d^{2}E$ for $t^{\prime}=-0.2$ in Fig. 8. In this case, $\chi_{c}>0$ for at least $x>0.06$. The phase separation area decreases for $t^{\prime}<0$.

Figure 8 indicates that $d^{2}E$ shows a depression near $x=0.12$ for $t^{\prime}<0$. This suggests the existence of strong charge fluctuation. We expect that this shows an instability toward some charge-ordered state such as the striped state.

In our optimized wave function, the instability toward the phase separation is limited to the range $x\equiv 1-n_{e}\leq 0.06$ for $t^{\prime}=0$ and the region of phase separation becomes small for negative $t^{\prime}$. We also mention that there is a possibility that the phase separation area will decrease as the wave function is optimized further by multiplying by operators $P_{G}$ and $\exp(-\lambda^{\prime}K)$.

5 Superconducting phase

Let us now discuss the SC state. We evaluate the SC condensation energy defined by

[TABLE]

where $\Delta_{SC}$ is the SC order parameter and $\Delta_{SC,{\rm opt}}$ is its optimized value. $\Delta E$ is the energy lowering due to the inclusion of the SC order parameter. We show the $U$-dependence of $\Delta_{SC,{\rm opt}}$ in Fig. 9 on a $10\times 10$ lattice for $N_{e}=88$ and $t^{\prime}=0$, where the upper curve is for the BCS-Gutzwiller function and the lower one is for the $\psi_{\lambda}$ function. Both curves in Fig. 9 have a maximum at $U/t\approx 12-14$. This shows that high-temperature superconductivity is possible in the strongly correlated region.

We discuss here the coexistence of SC and AF phases in the range of $0.06<x<0.09$. The wave function is written in the form

[TABLE]

The BCS wave function with both the SC order parameter $\Delta_{SC}$ and the AF order parameter $\Delta_{AF}$ is formulated by solving the Bogoliubov equation[42, 12]. We show the ground-state energy $E$ as a function of the electron number $N_{e}$ near $x\sim 0.08$ for $U/t=18$ in Fig. 10 for $\Delta_{SC}=0.005$, 0.01, 0.02, and 0.03, where the chemical potential $\mu$ in the BCS wave function is changed to adjust the number of electrons. The energy $E$ is lowered slightly by introducing $\Delta_{SC}$, by about $\Delta E/N\sim 0.005t$ per site at $\Delta_{SC}\sim 0.01t$. Here we used the parameters $g=0.0018$, $\lambda=0.055$, and $\Delta_{AF}=0.19t$.

The condensation energy per site as a function of the hole density (doping rate) $x$ is shown in Fig. 11 for $U/t=18$ and $t^{\prime}=0$ on a $10\times 10$ lattice. There is an instability toward the phase separation for $x\leq 0.06$. Thus, the AF state for $x\leq 0.06$ is an AF insulator. There is a coexistent metallic phase of superconductivity and antiferromagnetism when the doping rate is $0.06\leq x<0.09$. The SC condensation energy at $x=0.08$ is that for the coexistent state. In the range $0.09<x$, the pure $d$-wave SC state exists. We have also presented the result obtained by using the level-4 function with the AF order parameter in Fig. 11:

[TABLE]

where $g$, $g^{\prime}$, $\lambda$, and $\lambda^{\prime}$ are variational parameters. The condensation energy of the AF state for $\psi^{(4)}$ is less than that for $\psi_{\lambda}=\psi^{(2)}$.

With the inclusion of $t^{\prime}$, the phase separation region will decrease, and at the same time the area of the AF metallic state will increase. The phase diagram is dependent on $t^{\prime}$.

6 Summary

We have investigated the ground-state properties of the 2D Hubbard model by using the optimization variational Monte Carlo method. We used the exponential-type wave function given in the form $\exp(-\lambda S)$ with an appropriate operator $S$ and a variational parameter $\lambda$. With our wave function, the ground-state energy is lowered greatly and the energy expectation value is lower than that obtained by any other wave functions such as the Gutzwiller wave function and also several proposed wave functions with many variational parameters. The ground-state energy is lowered by the kinetic-energy gain originating from $\exp(-\lambda K)$.

The AF state is very stable near half-filling (with no carriers) in the 2D Hubbard model. The AF correlation is suppressed as the doping rate of holes increases. As the strength of the on-site Coulomb interaction $U$ increases, a crossover occurs between weakly correlated region and strongly correlated region. In the strongly correlated region, where $U$ is larger than $U_{c}$, which is of the order of the bandwidth, the AF correlation is suppressed. A decrease in the AF correlation indicates an increase in spin and charge fluctuation. This fluctuation is caused by an increase in kinetic energy and is likely induce electron pairing. We expect that this is an origin of high-temperature superconductivity.

We have shown the phase diagram in the plane of $\Delta E$ (condensation energy) and the hole doping rate $x$ for $t^{\prime}=0$. The value $t^{\prime}=0$ is most favorable for superconductivity, and for nonzero $t^{\prime}$ the AF area increases. In the underdoped region, where the doping rate is approximately $x<0.09$, AF order and superconductivity coexist.

We have also discussed the instability toward the phase separation in the low-doping region. The occurrence of this instability is dependent on the balance between the kinetic energy gain of holes and the electron interaction between the adjacent lattices such as the pairing interaction and the AF interaction. In the range $x\leq 0.06$, the state of phase separation is realized for $t^{\prime}=0$. In this region the AF state is an insulator, and becomes a metal when $x$ increases. The phase separation area decreases as $|t^{\prime}|$ increases. The region where AF and SC phases coexist becomes larger at the same time with increasing $-t^{\prime}$. We have also pointed out the depression in $\chi_{c}$ near $x=0.12$ for negative $t^{\prime}$ such as $t^{\prime}=-0.1\sim-0.2$.

This work was supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan (Grant No. 17K05559). The author expresses his sincere thanks to K. Yamaji and I. Hase for useful discussions. Part of the computations were supported by the Supercomputer Center of the Institute for Solid State Physics, the University of Tokyo.