Charge dynamics of correlated electrons formulated by variational description with inclusion of composite fermions

Charge dynamics of correlated electrons formulated by variational description
with inclusion of composite fermions

Kota Ido, Masatoshi Imada and Takahiro Misawa Institute for Solid State Physics, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan
Toyota Physical and Chemical Research Institute, Yokomichi, Nagakute, Aichi 480-1192, Japan
Research Institute for Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan

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 two-dimensional Hubbard models show that inclusion of composite-fermion 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 composite-fermion 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 single-particle descriptions manifested, for instance, by non-Fermi 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, pump-probe optical measurement, angle-resolved photoemission spectroscopy (ARPES), and resonant inelastic X-ray 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 X-ray 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 real-frequency () spectrum from the imaginary-time () quantities such as the Green’s function. It has been shown that this approach reasonably works for quantum-many body systemsSilver1990 (); Jarrell1996 (); Sandvik1998 (). The analytical continuation is, however, an ill-posed problem and sensitive to noises in the calculated imaginary-time data. In addition, the QMC method frequently suffers from the notorious negative-sign 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 many-body 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 real-frequency quantities are obtained from truncated density matrix for excitations. Another approach is the time-dependent DMRG (tDMRG) methodVidal2004 (); White2004 (); Daley2004 (). In the tDMRG, one first obtains the real-time dependence of physical quantities from the direct real-time 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 two-dimensional 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 time-dependent variational principleMclachlan1964 (); Haegeman2011 (), the imaginary-time or the real-time evolutions of the variational wave functions can be calculated. In principle, as is the case for the tDMRG approach, the real-frequency spectrum can be obtained by using the Fourier transformationsIdo2017 (). To perform the accurate Fourier transformation, however, the data of long real-time evolution is needed and the numerical cost becomes demanding. Therefore, a direct way of calculating the real-frequency 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 particle-hole 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 particle-hole excitations. One might think that the Li-Yang 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 Li-Yang 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 Li-Yang method with inclusion of only the simple particle-hole 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 Li-Yang method for dynamical quantities, which is based on the VMC method. Then we present our extension of composite-fermion 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 one-dimensional 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 two-dimensional 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 two-dimensional 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


Here, the matrix element of Hamiltonian and overlap matrices on the subspace are represented as


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:


Based on this approach, the dynamical spin structure factor of component of spin, can be computed. For example, is described as


where and represent the normalized ground-state 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


where and () is a creation operator for an electron with spin at position . Solution obtained only from this restricted Hilbert space is called the single-mode 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 particle-hole excitations is employed as the basis set of the restricted Hilbert space in the previous studies Li2010 (); DallaPiazza2015 (); Ferrari2018 (). This extension is given as


where represents the Gutzwiller factor to exclude the doubly occupied sites and is a mean-field solution for the ground state Ferrari2018 (). It was shown that the restricted basis set (9) can well describe the spinon continuum of in the one-dimensional Heisenberg model and in its extensionFerrari2018 ().

Advantages of this method are summarized as follows:

  1. No negative sign problems.

  2. Analytical continuation is not necessary.

  3. Excited states can be explicitly constructed from the ground-state 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 composite-fermion 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




In order to calculate in the VMC method, we approximate in the same manner as Eq. (II.1), namely


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 (equal-time) charge structure factor, namely

Next, we propose a way to restrict the basis set for . A naive candidate in analogy with would be


This approximation is the single-mode approximation of the density operator, , containing the long-range particle-hole 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 bare-fermion (BF) approach. At first glance, it looks possible to represent proper particle-hole 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 Mott-Hubbard 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


in addition to the bare electron operator


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 single-particle 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 composite-fermions as the composite-fermion (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


where and are the hopping integral and the on-site 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 Gutzwiller-JastrowGutzwiller1963 (); Jastrow1955 () correlation factors


The variational parameters , and are simultaneously optimized by using the stochastic reconfiguration methodSorella2001 (). We imposed the translational symmetry to the variational parameters for the one-dimensional case in order to reduce the numerical costs. This assumption works well for one-dimensional 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 doublon-holon 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

Figure 1: (Color online) Contour plot of the charge dynamical structure factor obtained by the VMC method in the two different approaches for the one-dimensional Hubbard model with , . The ED results are also shown for comparison. To enhance the visibility for difference between the results obtained by using two different variational Ansätze, the results are plotted in the logarithmic scale. The results for and are shown in panels (a-c) and (d-f), respectively. Panels (a) and (d) show the results obtained by using the ED method. In panels (b) and (e), we show the result obtained by the BF approach. The results by the CF approach are shown in panels (c) and (f). White dashed curves in panel (a) are dispersions characterizing the Bethe Ansatz solution: represents the upper (lower) edge energy of the two-holon [or two-doublon] continuum. represents the lower edge energy of the two-holon-two-spinon [or two-doublon-two-spinon] continuum. Notations of these modes are the same as those in the previous study Pereira2012 ().

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 particle-hole 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 two-holon (or two-doublon) continuum between and , and the other is the weak two-holon-two-spinon (or two-doublon-two-spinon) 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 two-holon-two-spinon continuum. This is reasonable because this continuum directly connects to particle-hole continuum in the limit of non-interacting 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 two-holon 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 two-holon continuum suggests that the BF approach does not describe the doublon-holon recombination process contributed from an excited doublon-holon 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 hole-doped 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 low-energy 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 low-energy 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 spin-charge 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.

Figure 2: (Color online) Charge dynamical structure factor in the one-dimensional Hubbard model for and at half filling for several choices of momentum . Black solid lines are the results obtained by using the ED method. Open squares and solid circles represent the VMC results obtained by using the basis set in Eq. (15) and Eq. (LABEL:new), respectively.
Figure 3: (Color online) Momentum-dependence of the absolute error of the charge dynamical structure factor in the one-dimensional Hubbard model for and at half filling. Open squares (solid circles) represents the results in the BF (CF) scheme.

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 many-body 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 lower-edge of the two-holon 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


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 hole-doped 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 low-energy 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.

Figure 4: (Color online) Contour plot of the charge dynamical structure factor in the one-dimensional Hubbard model for at half filling. The numerical method we used and the strength of the on-site interaction are described in the title of each panel. “ED” and “CF” represent the exact diagonalization method and the CF approach in the VMC method, respectively.
Figure 5: (Color online) Contour plot of the charge dynamical structure factor in the hole-doped Hubbard chain for and . The numerical method we used and the strength of the on-site interaction are described in the title of each panel. “ED” and “CF” denote the exact diagonalization method and the CF approach in the VMC method, respectively.

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 ().

Figure 6: (Color online) Contour plot of the charge dynamical structure factor obtained by using the CF approach in the Hubbard chain for and at half filling. The result is plotted in the logarithmic scale in the same way as Fig. 1. White dashed lines are dispersions characterizing the Bethe Ansatz solution as in Fig. 1 (a).

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 high-symmetry line for the system size . The boundary condition we used is the antiperiodic-periodic 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 nearest-neighbor doublon-holon 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 hole-doped 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.

Figure 7: (Color online) Contour plot of the charge dynamical structure factor in the two dimensional Hubbard model for and . is measured along a high-symmetry path through the Brillouin zone. The results for and are shown in panels (a-c) and (d-f), respectively. Panels (a) and (d) are obtained by using the exact diagonalization method. Panels (b) and (e) show the results obtained by using the BF approach defined in Eq. (15). Panels (c) and (f) show the results obtained by using the CF approach defined in Eq. (LABEL:new).

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


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 one-dimensional 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 higher-dimensional Mott insulators is an intriguing issue, which will be discussed in the next section in detail.

Figure 8: (Color online) Charge dynamical structure factor for (a) and (b) in the two-dimensional Hubbard model for and at half filling. Black solid lines are the results obtained by using the ED method. Open squares and solid circles represent the VMC results obtained by using the basis set in Eq. (15) and Eq. (LABEL:new), respectively.
Figure 9: (Color online) Interaction dependence of the mean absolute error in the Hubbard model at half filling. The dimension of the system and the VMC approach we used are shown in the legend.

For the hole-doped 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 low-energy 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.

Figure 10: (Color online) Charge dynamical structure factor for (a) and (b) in the two-dimensional Hubbard model for , and . Black solid lines are the results obtained by using the ED method. Open squares and solid circles represent the VMC results obtained by using the basis set in Eq. (15) and Eq. (LABEL:new), respectively.

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 spin-charge 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 ground-state 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 Jastrow-type 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 two-dimensional 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 low-energy 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 doublon-holon 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 long-ranged entanglement or correlation (from the nearest neighbour to several-distance neighbors) of the composite fermion such as with variational parameters and , defined as


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 many-body 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 one-particle spectral function . These extensions will be studied in the near future.

In summary, we examined how the VMC approach based on the Li-Yang method can be extended to describe the charge dynamics in the one- and two-dimensional Hubbard models. We found that the CF excitations are important for describing the charge dynamics in the Hubbard model. In the one-dimensional Hubbard model, we have shown that the CF approach well reproduces the results by the exact diagonalization at half filling and the hole-doped case as well. In the two-dimensional 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 two-dimensional 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 open-source 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 Grant-in-Aid 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 High-Performance 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




respectively. Therefore, we calculate up to sixth-order 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 node-difference among . See also Ref. Li2010, , where the reweighting technique is discussed in details.


  • (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. Sanchez-Palencia, 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. Tregenna-Piggott, 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).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description