Charge dynamics of correlated electrons formulated by variational description
with inclusion of composite fermions
Abstract
We propose a method to calculate the charge dynamical structure factors for the ground states of correlated electron systems based on the variational Monte Carlo method. Our benchmarks for the one and twodimensional Hubbard models show that inclusion of compositefermion excitations in the basis set greatly improves the accuracy, in reference to the exact charge dynamical structure factors for clusters. Together with examination for larger systems beyond tractable sizes by the exact diagonalization, our results indicate that the variational Monte Carlo method is a promising way for studies on the nature of charge dynamics in correlated materials such as the copper oxide superconductors, if the compositefermion excitations are properly included in the restricted Hilbert space of intermediate states in the linear response theory. We also discuss the importance of incorporating nonlocal composite fermion for more accurate description.
I Introduction
Strongly correlated electron systems are a platform to search the emergent properties of nature such as the breakdown of singleparticle descriptions manifested, for instance, by nonFermi liquid properties near Mott transitionsImada1998 (), quantum spin liquidsBalents2010 (); Zhou2017 (); Savary2017 (), and unconventional high superconductivityKeimer2015 (). To correctly understand the nature of the strongly correlated systems and control their physical properties, it is important to capture dominant part of elementary excitations emerging from correlations among electrons, which can be essentially different from what is expected from the noninteracting picture.
To experimentally measure those excitations, there exist several powerful techniques such as the neutron scattering, magnetic resonance, pumpprobe optical measurement, angleresolved photoemission spectroscopy (ARPES), and resonant inelastic Xray scattering (RIXS). To directly compare the experimental results obtained by the above techniques with theoretical estimates, it is necessary to develop accurate and efficient numerical methods to calculate dynamical physical quantities.
A seminal method for dynamical structure factors and dispersion of elementary bosonic excitations was formulated by Feynmann in his application to He, which succeeded in identifying rotons, in combination of neutron and Xray measurementFeynman1954 (). We follow basically the same spirit to calculate the dynamical structure factor by constructing excited states from the variational ground states satisfying the variational principle, now for strongly correlated fermionic systems. Thanks to the progress in computational power and methodology for dynamics since the Feyman’s work, the formulation is much more sophisticated than the original form of Feynmann as we describe below.
As one of the methods, analytical continuation of quantum Monte Carlo (QMC) data is often used for calculating dynamical structure factors. In this approach, by using the QMC method, one estimates the realfrequency () spectrum from the imaginarytime () quantities such as the Green’s function. It has been shown that this approach reasonably works for quantummany body systemsSilver1990 (); Jarrell1996 (); Sandvik1998 (). The analytical continuation is, however, an illposed problem and sensitive to noises in the calculated imaginarytime data. In addition, the QMC method frequently suffers from the notorious negativesign problem and its applicable range is limited. Therefore, an alternative approach without the analytical continuation and the sign problem is desired.
The density matrix renormalization group (DMRG) method is a powerful and accurate method for analyzing the quantum manybody systems in low dimensions without the sign problemsWhite1992 (); Schollwock2005 (). By using the DMRG, the dynamical quantities can be calculated efficiently in one dimensionHallberg1995 (); Ramasesha1997 (); Kuhner1999 (); Jeckelmann2000 (); Jeckelmann2002 (); Vidal2004 (); White2004 (); Daley2004 (). One of the widely used approaches is the dynamical DMRG (DDMRG) methodJeckelmann2000 (); Jeckelmann2002 (), where the realfrequency quantities are obtained from truncated density matrix for excitations. Another approach is the timedependent DMRG (tDMRG) methodVidal2004 (); White2004 (); Daley2004 (). In the tDMRG, one first obtains the realtime dependence of physical quantities from the direct realtime evolution. The real frequency quantities can be calculated after its Fourier transformation. However, the DMRG has difficulties in treating large systems in more than one dimension.
The variational Monte Carlo (VMC) method is another powerful method free from the sign problem becca2017quantum (). This method has been applied to analyze the ground states in a wide range of strongly correlated electron systems such as the twodimensional Hubbard model and its extensions Gros1987 (); Giamarchi1991 (); Sorella2001 (); Sorella2002 (); Yokoyama2013 (); Misawa2014a (); Tocchio2013 (); Zhao2017 (); Ido2018 (). Although the original applications of the VMC were limited to the ground states, a method of treating excited states and thermodynamic properties at nonzero temperatures were developed later by using approximated thermal pure statesTakai2016 () and transient states for nonequilibriumCarleo2012 (); Ido2015 (); Cevolani2015 (). In these approaches, based on the timedependent variational principleMclachlan1964 (); Haegeman2011 (), the imaginarytime or the realtime evolutions of the variational wave functions can be calculated. In principle, as is the case for the tDMRG approach, the realfrequency spectrum can be obtained by using the Fourier transformationsIdo2017 (). To perform the accurate Fourier transformation, however, the data of long realtime evolution is needed and the numerical cost becomes demanding. Therefore, a direct way of calculating the realfrequency properties within the VMC method is desired.
Recently, Li and Yang developed a VMC method to directly calculate the spin dynamical structure factorLi2010 (). The essence of their approach is to construct a relevant Hilbert subspace as the basis set by utilizing the Gutzwiller function applied to particlehole excitations generated from the approximated ground state. Within the restricted Hilbert subspace, dynamical properties can be directly calculated. Moreover, thanks to the Gutzwiller projection, their approach was shown to be capable of spinon dynamics in quantum spin chainsFerrari2018 () beyond simple particlehole excitations. One might think that the LiYang method for quantum spin dynamics could be straightforwardly extended to charge dynamics such as the dynamical charge structure factor
However, charge dynamics dominated by Mottness (barely itinerant electrons near the Mott insulator) together with interplay of spin fluctuations is a much more challenging subject. It involves competitions of charge, spin and superconducting order/fluctuations. If satisfactorily accurate methods are formulated, it has a wide range of applications in itinerant electron systems.
In this paper, we propose a nontrivial and efficient way of extension to the LiYang method to allow computation of the charge dynamical structure factors and test its accuracy by taking an example of the Hubbard model. Our benchmark calculations indicate that straightforward applications of the LiYang method with inclusion of only the simple particlehole excitations do not reproduce satisfactory charge dynamical structure factor even with Gutzwiller projection. By contrast, we show that the inclusion of the composite fermions into the basis set of the restricted subspace offers much better description of the charge dynamics. It is favorably compared with the exact diagonalization and the DMRG results.
The paper is organized as follows. In Sec. II, we first briefly review the essence of the LiYang method for dynamical quantities, which is based on the VMC method. Then we present our extension of compositefermion approach to calculate the charge dynamics. Section III presents our benchmarks for the charge dynamical structure factor for the Hubbard model. It is shown that introduction of the composite fermions drastically improves the charge dynamics in the Hubbard model. In the onedimensional Hubbard model, we find that our method quantitatively reproduces the results of the exact diagonalization. Our results for larger sized systems in one dimension are also consistent with the previous tDMRG resultPereira2012 (). Application to the twodimensional Hubbard model shows fairly good agreement with the exact results as well, representing charge dynamics of the undoped and doped Mott insulators including the incoherent part of dynamics. We, however, point out the possible further improvement by incorporating additional composite fermions not taken into account in this paper. Finally, in Sec. IV we summarize our paper and discuss future directions.
Ii Method
ii.1 Variational approach for excited states
In this subsection, we review a variational approach to calculate the dynamical spin structure factor . This method was introduced by Li and Yang to obtain in the  modelLi2010 (), and it has been recently applied to the one Ferrari2018 () and twodimensional Heisenberg model DallaPiazza2015 (); Ferrari2018a (). The idea of this method is to restrict the Hilbert space to a set of , which has momentum with an index which specifies type of excitations. In general, the basis set on the restricted Hilbert subspace is nonorthogonal and thus the generalized eigenvalue problem within this subspace can be written as
(1) 
Here, the matrix element of Hamiltonian and overlap matrices on the subspace are represented as
(2)  
(3) 
By solving this generalized eigenvalue problem defined in Eq. (1), we can obtain the th eigenvalue and the coefficients of its eigenvector . In other words, the excited state within the subspace is expressed in the following equation:
(4) 
Based on this approach, the dynamical spin structure factor of component of spin, can be computed. For example, is described as
(5)  
(6) 
where and represent the normalized groundstate wavefunction and its energy, respectively. Hereafter, we assume that the total momentum of the ground state is zero. We note that is a phenomenological smearing factor. By inserting the complete set of momentum in the restricted Hilbert space into Eq. (5), we obtain
To obtain accurate , choice of the basis set, namely how to pick up the restricted Hilbert space, is important. A simplest choice is
(8) 
where and () is a creation operator for an electron with spin at position . Solution obtained only from this restricted Hilbert space is called the singlemode approximation. In that case, one can easily compute the pole position because the dimension of the subspace at is only one. However, this approach can capture only an isolated dispersion, and is not able to represent continuum spectra if it exists as is expected in correlated systems. To overcome this limitation, the Gutzwiller wave function with particlehole excitations is employed as the basis set of the restricted Hilbert space in the previous studies Li2010 (); DallaPiazza2015 (); Ferrari2018 (). This extension is given as
(9) 
where represents the Gutzwiller factor to exclude the doubly occupied sites and is a meanfield solution for the ground state Ferrari2018 (). It was shown that the restricted basis set (9) can well describe the spinon continuum of in the onedimensional Heisenberg model and in its extensionFerrari2018 ().
Advantages of this method are summarized as follows:

No negative sign problems.

Analytical continuation is not necessary.

Excited states can be explicitly constructed from the groundstate wave function.
The third advantage can be rephrased as it is necessary to take proper component of the excited states into account to accurately calculate the dynamical quantities. As we show in the rest of this paper, incorporating the compositefermion excitations into the restricted Hilbert space is essential for accurate description of charge dynamical structure factors.
ii.2 Construction of excited states for charge dynamics
Following the idea in the previous studies, we propose a way to reasonably approximate the charge structure factor in the strongly correlated itinerant electron systems such as the Hubbard model. The definition of is
(10)  
(11) 
where
(12) 
In order to calculate in the VMC method, we approximate in the same manner as Eq. (II.1), namely
(13) 
We note that, in this scheme using the restricted Hilbert space, the sum rule is satisfied. Therefore, the integral of the approximated over is reduced to the static (equaltime) charge structure factor, namely
Next, we propose a way to restrict the basis set for . A naive candidate in analogy with would be
(15) 
This approximation is the singlemode approximation of the density operator, , containing the longrange particlehole excitation independent of spin degrees of freedom . Throughout this paper, the VMC approach where Eq. (15) is employed as the basis set of the restricted Hilbert space is called the barefermion (BF) approach. At first glance, it looks possible to represent proper particlehole continuum in as the case of . However, as we show later, this simple choice does not work well for describing the charge dynamics in the Mott insulating states. This is because holons and doublons represent essential charge states in the Mott insulating state of the Hubbard model, generating upper and lower Hubbard bands, which is not represented directly by bare electron operators and . The charge spectrum in the strong coupling limit is dominated by kinetics of excited holons and doublons in MottHubbard bands. To take this kinetics into account, we need to consider electron dynamics which depends on the occupation number of the sites in the process before and after the hopping. The importance of the separation of the electron kinetics has been already pointed out in the strong coupling approach based on the matrix expansion Gros1987 (); MacDonald1990 (). It is also related to the electron fractionalization, which is exactly proven in the strong coupling limit Zhu2013 (); Imada2019 ().
Thus, we need to distinguish these electron fractionalization effects beyond the bare electron to capture correct physical description of charge dynamics of correlated electrons. We therefore introduce composite fermion operators , defined as
(16)  
(17) 
in addition to the bare electron operator
(18) 
where denotes the opposite spin of . The composite fermions are known as the Hubbard operatorsHubbard1965 (); Mancini2004 (). Since creates a spinon (holon) and creates a doublon (spinon), they allow to differentiate dynamics of these particles. The importance of composite fermions to describe the nature of electronic structure in correlated materials has also been discussed in the structure of the singleparticle spectral functionYamaji2011 (); Yamaji2011a (); Zhu2013 (); Sakai2016a (); Sakai2016 (); Imada2019 (). We further note that can be expressed as and we need only two of these three because they are linearly dependent.
By using these operators, we here propose another candidate of the basis set as follows:
We call this VMC approach combined with the concept of compositefermions as the compositefermion (CF) approach. We evaluate the Hamiltonian matrix within this restricted Hilbert space and the overlap matrix on the target subspace, defined in Eqs. (2) and (3), by using the reweighting technique for efficient Monte Carlo samplingsLi2010 (). See also Appendix A for the detail of the reweighting technique.
Iii Results
iii.1 Model and setting
As benchmarks, we performed simulations of the dynamical structure factors in the Hubbard model. The Hubbard model is defined as
(20) 
where and are the hopping integral and the onsite Coulomb repulsion, respectively. In the following results, we set the energy unit to the hopping integral (). We employed the smearing factor for demonstration of . We note that there exists a trivial strong peak at for , whose integrated spectral weight is nothing but the square of the total number of particles, i.e, for . To eliminate this trivial peak and enhance visibility of the spectrum for other , we impose .
The trial wavefunction for the ground state we employed is a pair product wave functionTahara2008a () with GutzwillerJastrowGutzwiller1963 (); Jastrow1955 () correlation factors
(21)  
(22)  
(23)  
(24) 
The variational parameters , and are simultaneously optimized by using the stochastic reconfiguration methodSorella2001 (). We imposed the translational symmetry to the variational parameters for the onedimensional case in order to reduce the numerical costs. This assumption works well for onedimensional systems because the ground state should not spontaneously break the symmetry even in the thermodynamic limit. Note that our trial wavefunction can represent the Mott insulator owing to the doublonholon binding correlations in the Jastrow factorCapello2005 (). On the other hand, for the two dimensional case, we impose sublattice structure on to allow description of the antiferromagnetic states.
iii.2 One dimensional case
In this subsection, we show benchmarks of the charge dynamical structure factor in the one dimensional Hubbard model. We mainly study the system with linear size under the periodic boundary condition to directly compare with the exact diagonalization (ED) results. In the last part of this subsection, we show the result for as well to examine the applicable range of our method. At half filling, we employ a staggered particlehole transformation, and , in order to improve the accuracy of the ground stateShiba1972 (). We confirmed that the results do not change within the statistical error even if the transformation is not taken, although the statistical error is much higher.
First, to clarify whether composite fermions play an essential role to describe exact charge dynamics, we show the color contour plot of and dependence of obtained by two different variational Ansätze and compare with the ED results for in Fig. 1. At half filling, we draw three special modes obtained from the Bethe Ansatz equationPereira2012 (), and , in panel (a) by white dashed curves. In panel (a), there are two broad continuum at in the exact result. One is the strong twoholon (or twodoublon) continuum between and , and the other is the weak twoholontwospinon (or twodoublontwospinon) continuum between and for . Panel (b) shows the spectrum at half filling by using only the BF excitation in Eq. (15). We found that the BF approach is able to describe twoholontwospinon continuum. This is reasonable because this continuum directly connects to particlehole continuum in the limit of noninteracting system, which was mentioned in the previous DMRG calculation Pereira2012 (). However, the weight contributed from this continuum is larger than the exact one as we will detail in Fig. 2. In addition, this approach predicts a strong lower edge below the twoholon continuum. However, the ED results indicate that this region is dominated by a broad continuum especially for and no sharp edge is observed. The failure of capturing the broad twoholon continuum suggests that the BF approach does not describe the doublonholon recombination process contributed from an excited doublonholon pair. On the other hand, as we see in panel (c), the CF approach well reproduces the ED result.
Next, we show the results for the holedoped case in Fig. 1 (d)(f). Carrier doping eventually destroys the Mott gap, while the upper and lower Hubbard bands remain separated by the gap if the doping level is low. At low doping, the ground state becomes a correlated metal, where lowenergy excitation emerges near the Fermi level, but stays largely incoherent with broad feature particularly at large . In addition, although the incoherent excitation exists around as a remnant of the Mott gap at half filling, its weight is greatly reduced from the insulating case at higher doping level, which is transferred to the lowenergy part. The spectrum for finite doping calculated by using the BF approach is also shown in the panel (e). It underestimates the broadness at large . The reason would be that the hole carrier doping creates not only holons but also spinons due to the motion of the induced holons. Since the holon and spinon move separately because of the spincharge separation, we need to distinguish dynamics of these particles, which is beyond the representability of the single electron dynamics. On the other hand, the composite fermion scheme improves this broadness. The present results both for half filled and the hole doped cases indicate that the inclusion of the composite fermions are required to describe correct charge dynamics of correlated electrons.
The superiority of the CF description is more clearly shown in the following analysis. In Fig. 2 we show dependence of the charge structure factors for several momenta at half filling. At a small momentum () [shown in Fig. 2(a)], the spectrum obtained by the CF excitations well reproduces the ED result. The agreement of at the small momentum indicates that our CF approach has possibility to describe the correct optical conductivity as well because is tightly connected to the charge dynamical structure factor at , namely Stephan1996 (); Kim2004 (); Benthien2007 (). It is remarkable that the spinon excitations as a consequence of manybody effects in one dimension can essentially be captured by just two modes in the composite fermion scheme if the Gutzwiller and Jastrow factors are considered. The spectrum with the BF excitations, however, only reproduces the peak around and fails in representing other peaks. This contrast between the BF and CF is similar for larger [shown in Figs. 2(b),(c)]: The spectrum by BF excitations reproduces only one peak, which corresponds to the loweredge of the twoholon continuum . Furthermore, the amplitudes of the peak obtained by the BF excitations are larger than the exact results.
To analyze the difference between two approaches more quantitatively, Fig. 3 shows the absolute error of the charge dynamical structure factor for which is defined by
(25) 
Here is the charge dynamical structure factor obtained by using the exact diagonalization (VMC) method, and is the number of gird points on the line. We can confirm that by the CF approach is improved compared with that by the BF one.
Figures 4 and 5 show the interaction dependence of in the Hubbard chain at half filling and the holedoped case, respectively. In these figures, we only show by the ED and the CF approach because the accuracy of the BF approach turned out to be poor as we mentioned above. We see that in all the regions from weak to strong coupling, our VMC approach well reproduces the Mott gap scaled by , which is manifested in the global shift of the structure in the dependence of . For hole doped case, we find that the lowenergy excitation is not sensitive to the strength of , and the CF approach captures the dependence well, including the broad incoherent feature, suggesting the correct description of the holon dynamics. The relative weight of the upper Hubbard contribution at energy similar to the undoped case also shows good agreement with the ED results.
We also applied this method to a large size system that cannot be treated by the ED. The result is shown in Fig. 6. For large system size, the broad continuum around is clearly seen. We found that the strong lower edge and weak continuum below the edge for appear in the spectrum, which are consistent with the previous DMRG calculation Pereira2012 ().
iii.3 Two dimensional case
In this subsection, we present the results for the Hubbard model on the square lattice for . In Fig. 7, we show on the highsymmetry line for the system size . The boundary condition we used is the antiperiodicperiodic boundary condition to satisfy the closed shell condition at half filling.
We first examine the accuracy of our variational descriptions at half filling as shown in Figs. 7(a)(c). At half filling, the ED result shows that the strong but broad peak appears around due to the nearestneighbor doublonholon bindings. As is the case in one dimension, the BF approach does not describe charge dynamics even at half filling, i.e., the strong intensity at around surrounded by broad structure is not well reproduced. The CF approach clearly improves the description of this feature seen in the ED.
For holedoped case in two dimensions, as shown in Figs. 7(d)(f), we again found that the CF approach largely improves the charge dynamics compared with the BF approach. Especially, the broad continuum is well described by the CF approach in contrast to the BF approach. These results show that the CF approach works well for describing the charge dynamics even in two dimensions both at half filling and the doped case.
There is, however, still discrepancies from the exact results. A discrepancy is clearly seen in Figs. 8 (a) and (b), which shows dependence of for and , respectively. For , although only the single sharp peak appears around in the BF approach, this becomes broad structure in the CF approach and thus we can see the improvement by introducing the composite fermions. This improvement suggests that this broad structure in the ED result originates from the hybridization between the composite fermions defined in Eqs. (16) and (17). For , however, the improvement by introducing the composite fermions is insufficient: The broad spectrum around in the ED result is not captured in both of the VMC results. This suggests existence of other fermions, which are required to take into account to capture charge dynamics in the two dimensional Mott insulator as we discuss later.
This problem is universally seen irrespective of the strength of the Coulomb interaction . To see that, we plot the interaction dependence of the mean absolute error in Fig. 9, which is defined by
(26) 
where is the number of at which we measured . We see that the introduction of the composite fermions reduces for any strength of , However, for the two dimensional case obtained by the CF approach is still substantially larger than that for the onedimensional system. Our result indicates that the CF excitations may not be enough to describe the charge dynamics in the Mott insulators in higher dimensions. Identifying the key excitations to quantitatively describe the charge dynamics in the higherdimensional Mott insulators is an intriguing issue, which will be discussed in the next section in detail.
For the holedoped case, we see that the CF approach reproduces the overall charge dynamics. However, we notice that the intensities of by the CF approach are also stronger than those by the ED, which is also found in Fig. 10 where for and is plotted. The discrepancy for the incoherent part around , especially for , may share the same origin for the insufficient description of the charge dynamics in the Mott insulator at half filling as we have already shown. For , the spectral weight shows a lowenergy nearly flat dispersion in Fig. 7(f). On the other hand the exact result in Fig. 7(d) has much broader and damped feature. See also Fig. 10, which shows that the broadeness around for and for is underestimated and the intensity for and is overestimated. We discuss the possible origin of the discrepancy in the next section.
Iv Discussions and summary
The proposed VMC approach incorporating the composite fermion excitations provides us with accurate charge dynamical structure factors in one dimension and correctly reproduces signature of spincharge separation with broad incoherent continuum. It also shows a fairly good agreement with the ED results in two dimensions, with broad continuum and large intensity around in the incoherent part. However, we still have rooms for improvement in both the undoped case and the doped Mott insulator in two dimensions.
Here, we discuss the possible ways to improve the charge dynamics in two dimensions. There are two possible reasons for the insufficient broad continuum in the two dimensional system. One possible origin is the insufficient accuracy of the ground state, because the relative error of the groundstate energy tends to be larger than that in the chain Tahara2008a (); Kaneko2013 (). Recent proposed complementary methods for quantum lattice models such as the introductions of backflow correlations Tocchio2008 (); Tocchio2011 (); Ido2015 () and tensor network Chou2012 (); Sikora2015 (); Zhao2017 () could improve the Jastrowtype trial wavefunctions.
Second is the limitation of the assumed restricted basis set for the excitation in extracting the correct broad continuum. To enlarge the restricted subspace, we would need to introduce an efficient basis set for excited states in two dimensional systems. A simple candidate is another composite fermion which depends on intersite configurations. Such a composite fermion is required to consider charge dynamics in the antiferromagnetic background with spin fluctuations in the twodimensional systems, suggested in previous studiesAvella2003 (); Odashima2005 ().
The present construction of excited states is able to describe the upper and lower Hubbard excitations correctly because of the form Eqs.(16) and (17). However, in the doped system, ingap states are known to emerge, which also generates the pseudogap structure Ohta1992 (); Leung1992 (); Dagotto1992 (); Preuss1995 (); Preuss1997 (); Kyung2006 (); Sakai2009 (). Such lowenergy excitations near the Fermi level are not explicitly contained in the construction Eqs. (16) and (17) and likely to fail in capturing this emergent structure in the present form. It was proposed that such ingap states are generated by the coupling to excitonic excitations with weak binding energy in contrast to the doublonholon binding generating the large Mott gapImada2019 (). Such weakly bound excitons are moreover hybridizing with the spin singlet with the resonating valence bond nature Anderson1987 () and may give distinct ingap structure. This idea offers an alternative view to the coupling to the spin fluctuation in the antiferromagnetic background mentioned above. Both views suggest the importance to take into account longranged entanglement or correlation (from the nearest neighbour to severaldistance neighbors) of the composite fermion such as with variational parameters and , defined as
(27)  
(28)  
involving neighboring sites to as suggested in the context of the dark fermion in Ref.Imada2019, .
Both possibilities of improving the ground and excited states are desired to be examined toward more quantitative understanding of intriguing phenomena in quantum manybody systems. It is a fundamentally important issue for future studies.
We also note that the present approach would be also useful to obtain other dynamical physical quantities such as the optical conductivity and the oneparticle spectral function . These extensions will be studied in the near future.
In summary, we examined how the VMC approach based on the LiYang method can be extended to describe the charge dynamics in the one and twodimensional Hubbard models. We found that the CF excitations are important for describing the charge dynamics in the Hubbard model. In the onedimensional Hubbard model, we have shown that the CF approach well reproduces the results by the exact diagonalization at half filling and the holedoped case as well. In the twodimensional Hubbard model, although the CF approach largely improves the BF results, quantitative discrepancy from the exact results still exists. Our results indicate that the excitations complementing the present local CF excitations are necessary for quantitative description of the charge dynamics in the twodimensional Mott insulator. This is consistent with the emergence of the ingap states and the pseudogap, which are ignored in the present consideration of the local excitation. Spatially more extended and spin dependent composite fermions such as the dark fermion (or hidden fermion) Sakai2016 (); Sakai2018 (); Imada2019 () must be involved in the Hilbert space for the excitation. The role of coupling to weakly bound exciton is an intriguing future issue.
It is also an intriguing future study to apply our method to ab initio Hamiltonians for high superconductors Miyake2010 (); Hirayama2018 (); Nilsson2019 () and clarify how the enhanced dynamical charge fluctuations including the uniform static charge fluctuations affects high superconductivity Misawa2014a (); Misawa2014 (); Ohgoe2019 ().
Acknowledgments.— Our VMC code has been developed based on opensource software “mVMC”Misawa2017 (). To obtain the exact diagonalization results, we used “” packageKawamura2017 (). KI thanks Takahiro Ohgoe for continuous discussion. We acknowledge Maxime Charlebois and Youhei Yamaji for fruitful discussions and important suggestions. We thank the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo for the facilities. The authors are grateful to the support by a GrantinAid for Scientific Research (Nos. 16K17746, 16H06345, 19K03739, 19K14645) from Ministry of Education, Culture, Sports, Science and Technology, Japan. TM and KI were supported by Building of Consortia for the Development of Human Resources in Science and Technology from the MEXT of Japan. This work is financially supported by the MEXT HPCI Strategic Programs, and the Creation of New Functional Devices and HighPerformance Materials to Support Next Generation Industries (CDMSI). We also acknowledge the support provided by the RIKEN Advanced Institute for Computational Science under the HPCI System Research project (Grants No. hp170263, hp180170 and hp190145).
Appendix A Reweighting technique
In this appendix, we introduce the reweighting technique to efficiently evaluate the matrix element of Hamiltonian and overlap matrices, Eqs. (2) and (3) respectively, as proposed by Li and YangLi2010 (). By using the reweighting technique, we evaluate Eqs. (2) and (3) as
(30)  
(31)  
(32)  
(33) 
and
(34)  
(35) 
respectively. Therefore, we calculate up to sixthorder correlation functions to evaluate and at each sample. Since the weight is dependent on all bases we considered, we can reduce the statistical error caused by the nodedifference among . See also Ref. Li2010, , where the reweighting technique is discussed in details.
References
 (1) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
 (2) L. Balents, Nature 464, 199 (2010).
 (3) Y. Zhou, K. Kanoda, and T.K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
 (4) L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2017).
 (5) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
 (6) R. P. Feynman, Phys. Rev. 94, 262 (1954).
 (7) R. N. Silver, D. S. Sivia, and J. E. Gubernatis, Phys. Rev. B 41, 2380 (1990).
 (8) M. Jarrell and J. Gubernatis, Phys. Rep. 269, 133 (1996).
 (9) A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
 (10) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 (11) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
 (12) K. A. Hallberg, Phys. Rev. B 52, R9827 (1995).
 (13) S. Ramasesha, S. K. Pati, H. Krishnamurthy, Z. Shuai, and J. Brédas, Synth. Met. 85, 1019 (1997).
 (14) T. D. Kühner and S. R. White, Phys. Rev. B 60, 335 (1999).
 (15) E. Jeckelmann, F. Gebhard, and F. H. L. Essler, Phys. Rev. Lett. 85, 3910 (2000).
 (16) E. Jeckelmann, Phys. Rev. B 66, 045114 (2002).
 (17) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
 (18) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
 (19) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. Theory Exp. 2004, P04005 (2004).
 (20) F. Becca and S. Sorella, Quantum Monte Carlo approaches for correlated systems (Cambridge University Press, 2017).
 (21) C. Gros, R. Joynt, and T. M. Rice, Phys. Rev. B 36, 381 (1987).
 (22) T. Giamarchi and C. Lhuillier, Phys. Rev. B 43, 12943 (1991).
 (23) S. Sorella, Phys. Rev. B 64, 024512 (2001).
 (24) S. Sorella, G. B. Martins, F. Becca, C. Gazza, L. Capriotti, A. Parola, and E. Dagotto, Phys. Rev. Lett. 88, 117002 (2002).
 (25) H. Yokoyama, M. Ogata, Y. Tanaka, K. Kobayashi, and H. Tsuchiura, J. Phys. Soc. Jpn. 82, 014707 (2013).
 (26) T. Misawa and M. Imada, Phys. Rev. B 90, 115137 (2014).
 (27) L. F. Tocchio, H. Lee, H. O. Jeschke, R. Valentí, and C. Gros, Phys. Rev. B 87, 045111 (2013).
 (28) H.H. Zhao, K. Ido, S. Morita, and M. Imada, Phys. Rev. B 96, 085103 (2017).
 (29) K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 97, 045138 (2018).
 (30) K. Takai, K. Ido, T. Misawa, Y. Yamaji, and M. Imada, J. Phys. Soc. Jpn. 85, 034601 (2016).
 (31) G. Carleo, F. Becca, M. Schiró, and M. Fabrizio, Sci. Rep. 2, 243 (2012).
 (32) K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 92, 245106 (2015).
 (33) L. Cevolani, G. Carleo, and L. SanchezPalencia, Phys. Rev. A 92, 041603(R) (2015).
 (34) A. McLachlan, Mol. Phys. 8, 39 (1964).
 (35) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
 (36) K. Ido, T. Ohgoe, and M. Imada, Sci. Adv. 3, e1700718 (2017).
 (37) T. Li and F. Yang, Phys. Rev. B 81, 214509 (2010).
 (38) F. Ferrari, A. Parola, S. Sorella, and F. Becca, Phys. Rev. B 97, 235103 (2018).
 (39) R. G. Pereira, K. Penc, S. R. White, P. D. Sacramento, and J. M. P. Carmelo, Phys. Rev. B 85, 165132 (2012).
 (40) B. Dalla Piazza, M. Mourigal, N. B. Christensen, G. J. Nilsen, P. TregennaPiggott, T. G. Perring, M. Enderle, D. F. McMorrow, D. A. Ivanov, and H. M. Rønnow, Nat. Phys. 11, 62 (2015).
 (41) F. Ferrari and F. Becca, Phys. Rev. B 98, 100405(R) (2018).
 (42) A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 41, 2565 (1990).
 (43) L. Zhu and J.X. Zhu, Phys. Rev. B 87, 085120 (2013).
 (44) M. Imada and T. J. Suzuki, J. Phys. Soc. Jpn. 88, 024701 (2019).
 (45) J. Hubbard, Proc. Roy. Soc. A 285, 542 (1965).
 (46) F. Mancini and A. Avella, Adv. Phys. 53, 232 (2004).
 (47) Y. Yamaji and M. Imada, Phys. Rev. B 83, 214522 (2011).
 (48) Y. Yamaji and M. Imada, Phys. Rev. Lett. 106, 016404 (2011).
 (49) S. Sakai, M. Civelli, and M. Imada, Phys. Rev. Lett. 116, 057003 (2016).
 (50) S. Sakai, M. Civelli, and M. Imada, Phys. Rev. B 94, 115130 (2016).
 (51) D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77, 114701 (2008).
 (52) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
 (53) R. Jastrow, Phys. Rev. 98, 1479 (1955).
 (54) M. Capello, F. Becca, M. Fabrizio, S. Sorella, and E. Tosatti, Phys. Rev. Lett. 94, 026406 (2005).
 (55) H. Shiba, Prog. Theor. Phys. 48, 2171 (1972).
 (56) W. Stephan and K. Penc, Phys. Rev. B 54, R17269 (1996).
 (57) Y.J. Kim, J. P. Hill, H. Benthien, F. H. L. Essler, E. Jeckelmann, H. S. Choi, T. W. Noh, N. Motoyama, K. M. Kojima, et al., Phys. Rev. Lett. 92, 137402 (2004).
 (58) H. Benthien and E. Jeckelmann, Phys. Rev. B 75, 205128 (2007).
 (59) R. Kaneko, S. Morita, and M. Imada, J. Phys. Conf. Ser. 454, 012046 (2013).
 (60) L. F. Tocchio, F. Becca, A. Parola, and S. Sorella, Phys. Rev. B 78, 041101(R) (2008).
 (61) L. F. Tocchio, F. Becca, and C. Gros, Phys. Rev. B 83, 195138 (2011).
 (62) C.P. Chou, F. Pollmann, and T.K. Lee, Phys. Rev. B 86, 041105(R) (2012).
 (63) O. Sikora, H.W. Chang, C.P. Chou, F. Pollmann, and Y.J. Kao, Phys. Rev. B 91, 165113 (2015).
 (64) A. Avella, F. Mancini, and T. Saikawa, Eur. Phys. J. B 36, 445 (2003).
 (65) S. Odashima, A. Avella, and F. Mancini, Phys. Rev. B 72, 205121 (2005).
 (66) Y. Ohta, K. Tsutsui, W. Koshibae, T. Shimozato, and S. Maekawa, Phys. Rev. B 46, 14022 (1992).
 (67) P. W. Leung, Z. Liu, E. Manousakis, M. A. Novotny, and P. E. Oppenheimer, Phys. Rev. B 46, 11779 (1992).
 (68) E. Dagotto, F. Ortolani, and D. Scalapino, Phys. Rev. B 46, 3183 (1992).
 (69) R. Preuss, W. Hanke, and W. von der Linden, Phys. Rev. Lett. 75, 1344 (1995).
 (70) R. Preuss, W. Hanke, C. Gröber, and H. G. Evertz, Phys. Rev. Lett. 79, 1122 (1997).
 (71) B. Kyung, S. S. Kancharla, D. Sénéchal, A.M. S. Tremblay, M. Civelli, and G. Kotliar, Phys. Rev. B 73, 165114 (2006).
 (72) S. Sakai, Y. Motome, and M. Imada, Phys. Rev. Lett. 102, 056404 (2009).
 (73) P. W. Anderson, Science 235, 1196 (1987).
 (74) S. Sakai, M. Civelli, and M. Imada, Phys. Rev. B 98, 195109 (2018).
 (75) T. Miyake, K. Nakamura, R. Arita, and M. Imada, J. Phys. Soc. Jpn. 79, 044705 (2010).
 (76) M. Hirayama, Y. Yamaji, T. Misawa, and M. Imada, Phys. Rev. B 98, 134501 (2018).
 (77) F. Nilsson, K. Karlsson, and F. Aryasetiawan, Phys. Rev. B 99, 075135 (2019).
 (78) T. Misawa and M. Imada, Nat. Commun. 5, 5738 (2014).
 (79) T. Ohgoe, M. Hirayama, T. Misawa, K. Ido, Y. Yamaji, and M. Imada, arXiv:1902.00122 .
 (80) T. Misawa, S. Morita, K. Yoshimi, M. Kawamura, Y. Motoyama, K. Ido, T. Ohgoe, M. Imada, and T. Kato, Comput. Phys. Commun. 235, 447 (2019).
 (81) M. Kawamura, K. Yoshimi, T. Misawa, Y. Yamaji, S. Todo, and N. Kawashima, Comput. Phys. Commun. 217, 180 (2017).