# Variational Monte Carlo Study of Electron Differentiation around Mott Transition

\abstWe study ground-state properties of the two-dimensional Hubbard model at half filling by improving variational Monte Carlo method and by implementing quantum-number projection and multi-variable optimization. The improved variational wave function enables a highly accurate description of the Mott transition and strong fluctuations in metals. We clarify how anomalous metals appear near the first-order Mott transition. The double occupancy stays nearly constant as a function of the on-site Coulomb interaction in the metallic phase near the Mott transition in agreement with the previous unbiased results. This unconventional metal at half filling is stabilized by a formation of “electron-like pockets” coexisting with an arc structure, which leads to a prominent differentiation of electrons in momentum space. An abrupt collapse of the “pocket” and “arc” drives the first-order Mott transition. \kwordvariational Monte Carlo method, strongly correlated electron systems, Hubbard model, Mott transition, Fermi arc When electrons become localized by the strong Coulomb interaction, the Mott insulator appears after a metal-insulator transition called the Mott transition [1, 2]. The Mott transition is found in many materials and systems, such as transition metal oxides [2], organic conductors [3], and He layered systems [4, 5].

Filling-control Mott transitions at zero temperature have been studied by the auxiliary-field quantum Monte Carlo (AFQMC) method in the Hubbard model on a square lattice [6, 7]. The transition shows a continuous character with divergences of the compressibility and the antiferromagnetic correlation length. Bandwidth-control Mott transitions have been studied by the path-integral renormalization group (PIRG) method [8, 9, 10]. The results have shown the first-order as well as continuous Mott transitions. The Mott transition has also been studied by the dynamical mean-field theory (DMFT) [11]. The DMFT has shown that the Mott transition at the critical end point at nonzero temperature is consistent with the Ising universality class [12]. Extensions of the DMFT to include spatial correlations [13, 14, 15], such as the cellular dynamical mean-field theory (CDMFT), have also suggested the existence of both types of the first-order and continuous transitions. The Mott transition at finite temperatures has been studied by a phenomenological effective theory [16] and a mean-field theory as well [17, 18]. They have found unconventional criticalities at a marginal quantum critical point (MQCP), which arises from an interplay of topological transition and symmetry breaking. Experimentally measured critical exponents of an organic conductor, -(BEDT-TTF)Cu[N(CN)]Cl are in agreement with those of the MQCP [19].

While crucial properties of the Mott transition, such as criticalities, have been elucidated by recent extensive studies, underlying electronic structure of metals near the Mott insulators has not fully been understood. This is an important issue to be solved, because competing orders and quantum phases near the Mott transition must be consequences of the underlying unconventional electronic structure. Filling-control Mott transitions generating the high-temperature superconductivity in the copper oxides, and bandwidth-control Mott transitions generating the superconductivity and the quantum spin liquid in the organic conductors [20] are such consequences.

Recently, electron differentiations in momentum space around the Mott transition are extensively studied both in experimental [21] and theoretical [22, 23] studies as a key to understand the unconventional electronic structure. The angle resolved photoemission spectroscopy has revealed a truncation of the Fermi surface called Fermi arc, in the underdoped region of the high-temperature copper oxide superconductors (HTSC) [21]. This truncation is a typical differentiation where electrons at different points of the Fermi surface start showing distinct behavior. The origin of the Fermi arc and the pseudogap in the HTSC has been extensively studied by using the extension of DMFT, such as the CDMFT [24, 25]. In order to capture the differentiation, an accurate treatment of spatial correlations and high resolution in momentum space are crucially important, while extensions of CDMFT to larger spatial sizes are difficult. For this purpose, the variational Monte Carlo (VMC) method [26] offers an alternative advantageous approach, where large system sizes are tractable even with strong interactions and geometrical frustration effects. It offers high-resolution results in momentum space. However, the bias inherently and inevitably contained in the assumed variational form of the wave functions is a severe limitation in the VMC method. Reducing the biases is crucially important in studies with variational wave functions.

In this paper, we study the electron differentiation of the Mott transition by improving the VMC method. We apply our improved variational wave function [27] to the two-dimensional Hubbard model with geometrical frustration effects. We find that highly fluctuating metals with large amplitude of double occupancy of electrons is retained near the Mott transition. This unconventional metal at half filling is stabilized by the electron differentiation through a formation of “electron-like pockets” together with an arc structure.

The Hubbard model on a square lattice is defined as

(1) |

where () is the creation (annihilation) operator on the -th site with spin and is the number operator. The transfer integral is taken as for the nearest-neighbor sites and for the next nearest-neighbor sites. We take sites with the boundary condition periodic in the direction and antiperiodic in the direction (periodic-antiperiodic boundary condition). Throughout this paper, we consider the half-filled case at .

The variational wave function used in this study has the following form [27]:

(2) |

where , , and are the Gutzwiller factor [28], the doublon-holon correlation factor [29, 30], and the Jastrow factor [31], respectively. These factors are defined as

(3) | ||||

(4) | ||||

(5) |

where , , and are variational parameters. We assume that depends on the displacement . Here, is a many-body operator which is diagonal in the real space representations. When a doublon (holon) exists at the -th site and holons (doublons) surround at the -th nearest neighbor, returns . In other cases, returns . The spin quantum-number projection restores the spin-rotational symmetry and generates a state with total spin [32, 33]. The one-body part is the generalized pairing wave function defined as

(6) |

with the conditions and . Here, BZ and AFBZ denote the Brillouin zone and the folded antiferromagnetic (AF) Brillouin zone, respectively. The pair orbitals and are variational parameters. Since the PIRG method shows that the ordering vector of the AF insulating phase in the interval, , is [8, 10], the vector in eq. (6) is chosen as in this study. All the variational parameters (, , , , and ) are simultaneously optimized by using the stochastic reconfiguration method [34]. The variational wave function in eq. (2) allows to describe various states such as paramagnetic metals (PM), antiferromagnetically ordered states, and superconducting states with any wavenumber (spatial) dependence of gap functions within a single functional form. Moreover, enables efficient treatment of quantum fluctuations with long-ranged as well as short-ranged correlations. For example, the peak value of the spin structure factor in the doped Hubbard model shows quantitative agreement with the results obtained from the unbiased AFQMC method. Detailed discussions and extensive benchmarks are reported elsewhere [27].

In order to capture the metal-insulator and magnetic phase transitions, we calculate the double occupancy , spin structure factor , and charge gap for -site systems, where is the spin operator and the chemical potential is given as . Here, is the variational energy with the number of up (down) spin (). The double occupancy in the thermodynamic limit is estimated by fitting the finite-size data in the form , because the PIRG results for the frustrated Hubbard model have succeeded in the extrapolation by this form and estimating the critical point for the MIT [8, 9]. The staggered magnetization in the thermodynamic limit is also estimated by fitting the scaling form with by following the scaling of the spin wave theory if the long-range order is present [35]. We also calculate the local moment defined as . For the charge gap in the thermodynamic limit , we use the scaling function as in the AF Hartree-Fock gap equation [6]. The extrapolation to the thermodynamic limit is performed by using several choices of system sizes up to .

Figure 1 shows the dependence of , , , and . The first-order metal-insulator transition takes place at . The magnetic phase transition takes place at the same critical point within our resolution. The nonmagnetic insulating phase is not clearly identified in our variational results. Although gradually increases as a function of , the growth is strongly suppressed in the metallic phase near the metal-insulator transition. In the metallic region, at rapidly grows when approaches , although the true long-ranged order is not achieved within metals. This growth is, however, to a large extent, ascribed to the growth in the longer ranged part in the real space correlation, while the shorter ranged part of antiferromagnetic correlations does not grow equally. In this overall tendency, the local moment , that is the shortest-ranged part, shows nearly flat dependence on . This seems to optimize and reconcile in gaining the kinetic energy by keeping the electron-pocket-like and arc-like structure of the Fermi surface as we describe below. In other words, the flat dependence of with growing means a compensating reduction of at . This reduction suppresses the electron scattering by the spin fluctuations and keeps the coherence of the carrier, leading to the gain in the kinetic energy.

Figure 1 shows remarkable quantitative agreement with the unbiased results obtained from the PIRG method [8]. Although the value in the PIRG data in ref. 8 is not the same as that of the present results, the values at are rather precisely interpolated from the results at and given by the PIRG method. For example, estimated by interpolation of the PIRG results suggests , which may be compared with the present estimate . The absence or only tiny interval of the nonmagnetic insulating region at is also consistent each other. Furthermore, near in the metallic phase stays nearly constant around over an interval contrary to a naive expectation. This flat dependence of also remarkably agrees with the PIRG results. In addition, Fig. 1 shows convex growth of the squared staggered-magnetization from for . These features have been observed by the PIRG method at as well. In our results, the squared staggered-magnetization at and is . In the PIRG results, the value at and is . When we consider the difference of , this is a remarkable quantitative agreement. Our results suggest that our variational wave function enables quantitatively accurate descriptions in the ground state.

On the other hand, in the previous VMC calculation for this model, the critical value of the Mott transition is much larger () [36]. In addition, the double occupancy in ref. 36 linearly decays to the transition point as a function of with a substantial slope (). The variational wave functions employed in the literature include many-body correlations only in much restricted forms, such as the short-ranged doublon-holon factor. Such a restricted form does not sufficiently take into account fluctuations, which are strongly enhanced around the Mott transition. By introducing a large number of variational parameters that scales linearly with the system size, the Gutzwiller-Jastrow factor as well as the one-body part allow much more accurate treatment of fluctuations . DMFT results also show large even for and fail in capturing the plateau of at [11]. Although the CDMFT takes into account spatial correlations to some extent, the plateau of around has not been captured yet in the CDMFT studies [37, 38]. The exact diagonalization study [39] does not capture this behavior either because of the limitation of the system size. Sufficiently large system size over the correlation length of the fluctuations is important for reproducing the sufficient coherence and the plateau of .

In order to analyze the electron differentiation in this region,
we calculate the momentum distribution
.
We assume that has fourfold symmetric structure.
Since we calculate under the periodic-antiperiodic boundary condition,
this assumption allows twice as many -points in the BZ for the system:

(7) |

for . The VMC results for are shown in Fig. 2. The momentum distribution shows a characteristic behavior near the Mott transition. In the metallic phase, the amplitude of around points () is kept large, generating an “electron pocket-like” structure just before the Mott transition as is seen in Figs. 2(a) and (b).

Figures 3(a)-(d) are contour plots of . We can estimate the Fermi surface and the coherence of electrons on the Fermi surface by the amplitude of . The contour plots remarkably show the “arc-like” structure in the metallic phase close to the Mott transition (Fig. 3(c)). Although our system size is not sufficient, the position of the Fermi surface looks rigid and very similar to that of the non-interacting system, contrary to the “deformation to the nesting” obtained in the renormalization-group method in the weak-coupling regime [40]. The rigidity enhances the coexistence of “pocket” and “arc.” The “arc” around coexisting with the “electron pocket-like” structure around and discussed before is retained near the Mott transition. The coexistence implies that, although the long-range magnetic order is absent, the gap structure of the charge excitation already gets formed even in the metallic phase as a precursor of the Mott gap. As in semimetals, the preformed Mott gap generates electron and hole pockets [18]. The arc structure emerges as a part of the hole pockets around where outer half of the pocket is lost presumably because of a strong damping. The coexistence of the “electron pocket” and “arc” directly causes largely retained and suppressed in Fig. 1. This “semimetallic pocket” structure creates electrons and holes in the upper and lower Hubbard bands, respectively. As a result, excitons (or doublon-holon pairs) are generated in the background of the Mott insulator, where is largely retained. The first-order Mott transition appears as a sudden collapse of the “semimetallic pocket” structure. The similar abrupt change of is also seen in the PIRG results [8]. The Mott transition emerging as shrinkage of electron-hole (or doublon-holon) pockets supports the topological character of the transition proposed in ref. 18.

In summary, we have studied the electron differentiation in momentum space around the Mott transition by using the improved VMC method. Our variational results show quantitative agreement with the unbiased method. Especially, our variational wave function enables to treat short-ranged as well as long-ranged spin/charge fluctuations, which are strongly enhanced near the Mott transition. As a result, we have captured the flat dependence of double occupancy and abrupt change of momentum distribution at the first-order Mott transition. The momentum distribution shows “electron pocket-like” structure around . The arc structure appears around in in the metallic phase near the Mott transition. The abrupt collapse of these parts drives the first-order Mott transition. The coexistence of “electron pocket-like” structure and “hole-like arcs” with the retained plateau of the double occupancy is the key aspects of the electron differentiation in momentum space. Clarifying the relation between the above differentiation and instability toward superconductivity is one of the most important issues left for future studies.

One of the authors (D.T.) thanks S. Watanabe, T. Misawa, and Y. Yamaji for useful discussions. This work was supported by Grants-in-Aid for Scientific Research on Priority Areas under the grant numbers 17071003, 16076212, and 17064004 from the Ministry of Education, Culture, Sports, Science and Technology. A part of our computation has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

## References

- [1] N. F. Mott and R. Peierls: Proc. Soc. London. 49 (1937) 72.
- [2] For a review see M. Imada, A. Fujimori, and Y. Tokura: Rev. Mod. Phys. 70 (1998) 1039.
- [3] K. Kanoda: J. Phys. Soc. Jpn. 75 (2006) 051007.
- [4] K. Ishida, M. Morishita, K. Yawata, and H. Fukuyama: Phys. Rev. Lett. 79 (1997) 3451.
- [5] R. Masutomi, Y. Karaki, and H. Ishimoto: Phys. Rev. Lett. 92 (2004) 025301.
- [6] N. Furukawa and M. Imada: J. Phys. Soc. Jpn. 61 (1992) 3331.
- [7] N. Furukawa and M. Imada: J. Phys. Soc. Jpn. 62 (1993) 2557.
- [8] T. Kashima and M. Imada: J. Phys. Soc. Jpn. 70 (2001) 3052.
- [9] H. Morita, S. Watanabe, and M. Imada: J. Phys. Soc. Jpn. 71 (2002) 2109.
- [10] T. Mizusaki and M. Imada: Phys. Rev. B 74 (2006) 014421.
- [11] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg: Rev. Mod. Phys. 68 (1996) 13.
- [12] G. Kotliar, E. Lange, and M. J. Rozenberg: Phys. Rev. Lett. 84 (2000) 5180.
- [13] S. Onoda and M. Imada: Phys. Rev. B 67 (2003) 161102.
- [14] O. Parcollet, G. Biroli, and G. Kotliar: Phys. Rev. Lett. 92 (2004) 226402.
- [15] M. Aichhorn, E. Arrigoni, M. Potthoff, and W. Hanke: Phys. Rev. B 76 (2007) 224509.
- [16] M. Imada: Phys. Rev. B 72 (2005) 075113.
- [17] T. Misawa, Y. Yamaji, and M. Imada: J. Phys. Soc. Jpn. 75 (2006) 083705.
- [18] T. Misawa and M. Imada: Phys. Rev. B 75 (2007) 115121.
- [19] F. Kagawa, K. Miyagawa, and K. Kanoda: Nature 436 (2005) 534.
- [20] Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito: Phys. Rev. Lett. 91 (2003) 107001.
- [21] A. Damascelli, Z. Hussain, and Z.-X. Shen: Rev. Mod. Phys. 75 (2003) 473.
- [22] S. Onoda and M. Imada: J. Phys. Chem. Solids 62 (2001) 47.
- [23] S. Onoda and M. Imada: J. Phys. Chem. Solids 63 (2002) 2225.
- [24] T. D. Stanescu and G. Kotliar: Phys. Rev. B 74 (2006) 125110.
- [25] Y. Z. Zhang and M. Imada: Phys. Rev. B 76 (2007) 045108.
- [26] D. Ceperley, G. V. Chester, and M. H. Kalos: Phys. Rev. B 16 (1977) 3081.
- [27] D. Tahara and M. Imada: arXiv:0805.4457.
- [28] M. C. Gutzwiller: Phys. Rev. Lett. 10 (1963) 159.
- [29] T. A. Kaplan, P. Horsch, and P. Fulde: Phys. Rev. Lett. 49 (1982) 889.
- [30] H. Yokoyama and H. Shiba: J. Phys. Soc. Jpn. 59 (1990) 3669.
- [31] R. Jastrow: Phys. Rev. 98 (1955) 1479.
- [32] P. Ring and P. Schuck: The Nuclear Many-Body Problem, (Springer-Verlag, New York, Heidelberg, Berlin, 1980).
- [33] T. Mizusaki and M. Imada: Phys. Rev. B 69 (2004) 125110.
- [34] S. Sorella: Phys. Rev. B 64 (2001) 024512.
- [35] D. A. Huse: Phys. Rev. B 37 (1988) 2380.
- [36] H. Yokoyama, M. Ogata, and Y. Tanaka: J. Phys. Soc. Jpn. 75 (2006) 114706.
- [37] B. Kyung and A.-M. S. Tremblay: Phys. Rev. Lett. 97 (2006) 046402.
- [38] T. Ohashi, T. Momoi, H. Tsunetsugu, and N. Kawakami: Phys. Rev. Lett. 100 (2008) 076402.
- [39] T. Koretsune, Y. Motome, and A. Furusaki: J. Phys. Soc. Jpn. 76 (2007) 074719.
- [40] T. Ogawa, H. Maebashi, H. Kohno, and K. Miyake: Physica B 312-313 (2002) 525.