Single-particle spectral function formulated and calculated by variational Monte Carlo method with application to $d$-wave superconducting state

A method to calculate the one-body Green's function for ground states of correlated electron materials is formulated by extending the variational Monte Carlo method. We benchmark against the exact diagonalization (ED) for the one- and two-dimensional Hubbard models of 16 site lattices, which proves high accuracy of the method. The application of the method to larger-sized Hubbard model on the square lattice correctly reproduces the Mott insulating behavior at half filling and gap structures of $d$-wave superconducting state of the hole doped Hubbard model in the ground state optimized by enforcing the charge uniformity, evidencing a wide applicability to strongly correlated electron systems. From the obtained $d$-wave superconducting gap of the charge uniform state, we find that the gap amplitude at the antinodal point is several times larger than the experimental value, when we employ a realistic parameter as a model of the cuprate superconductors. The effective attractive interaction of carriers in the $d$-wave superconducting state inferred for an optimized state of the Hubbard model is as large as the order of the nearest-neighbor transfer, which is far beyond the former expectation in the cuprates. We discuss the nature of the superconducting state of the Hubbard model in terms of the overestimate of the gap and the attractive interaction in comparison to the cuprates.

A method to calculate the one-body Green's function for ground states of correlated electron materials is formulated by extending the variational Monte Carlo method. We benchmark against the exact diagonalization (ED) for the one-and two-dimensional Hubbard models of 16 site lattices, which proves high accuracy of the method. The application of the method to larger-sized Hubbard model on the square lattice correctly reproduces the Mott insulating behavior at half filling and gap structures of d-wave superconducting state of the hole doped Hubbard model in the ground state optimized by enforcing the charge uniformity, evidencing a wide applicability to strongly correlated electron systems. From the obtained d-wave superconducting gap of the charge uniform state, we find that the gap amplitude at the antinodal point is several times larger than the experimental value, when we employ a realistic parameter as a model of the cuprate superconductors. The effective attractive interaction of carriers in the d-wave superconducting state inferred for an optimized state of the Hubbard model is as large as the order of the nearest-neighbor transfer, which is far beyond the former expectation in the cuprates. We discuss the nature of the superconducting state of the Hubbard model in terms of the overestimate of the gap and the attractive interaction in comparison to the cuprates.

I. INTRODUCTION
Dynamical properties often provide us with crucial insights into open issues of strongly correlated electron systems. In particular, momentum and energy resolved single-particle spectral function A(k, ω), which is the imaginary part of the Green's function G(k, ω) with the momentum k and the frequency ω gives us understanding how an electron moves in an environment of other mutually interacting electrons and provides us with properties of the excited states, which in turn reveals the equilibrium properties as well.
In an example of the copper oxide high-T c superconductors, A(k, ω) has been extensively studied by the angle-resolved photoemission spectroscopy (ARPES), which has greatly contributed to elucidate the properties of superconducting as well as anomalous normal metallic properties including the pseudogap, Fermi arc and d-wave superconducting gap structure itself [1,2].
Numerical methods to clarify the dynamics of the strongly correlated electron systems have been hampered by various difficulties such as fermion sign problem [3][4][5][6] in quantum Monte Carlo methods, and intrinsic quantum entanglement of electrons at long distances. Nevertheless, linear response quantities such as the spin and charge dynamical structure factors, S(k, ω), N (k, ω), respectively, defined below have been studied by limited methods such as the exact diagonalization (ED) [7,8] and time-dependent density matrix renormalization group [9,10]. However, these methods have their own limitations, namely amenable only in small system sizes and in one-dimensional lattice structure, respectively.
In addition to the exact diagonalization, single-particle Green's function G(k, ω) has been studied by the cluster extension of the dynamical mean-field theory (cDMFT) [11][12][13], but the allowed momentum resolution is severely limited by the cluster size at most in the order of 10. The cDMFT method is generally combined with a periodization procedure in order to restore translation invariance and provides us with the data at interpolated momenta k [14][15][16][17]. However, we need to be cautious about these periodizations and the results should be regarded as estimators, because the momentum resolution is limited by the cluster size anyhow. This is true even for inhomogenous extension of cDMFT [18,19], where large supercluster still retains the self-energy modulation of the original smallest cluster [20]. The need to study bigger cluster remains.
Recent formulation of time-dependent variational Monte Carlo method based on the variational principle opened a way to study the long-time dynamics [21,22], but it has not been extensively applied yet to interacting fermion systems except for few examples [23]. Meanwhile, methods of calculating the spin and charge dynamical structure factors utilizing the variational wave functions for ground and excited states have been proposed recently [24][25][26][27][28]. Some attempts to calculate excitation spectrum on larger cluster of the t-J model [29,30].
Here, we formulate a method of calculating the Green's function G(k, ω) and the spectral function A(k, ω) = − 1 π ImG(k, ω) and show its accuracy by comparing with the exact results. It reproduces the feature of the spincharge separation and excitation continuum in the onedimensional Hubbard model.
The method is applied to the Hubbard model on the square lattice as well. In the optimized d-wave superconducting solution, though it is an excited state in the competition with the stripe states [31][32][33], the d-wave symmetry of the gap structure is correctly reproduced in the spectral function in this charge-uniform lowest energy state. However, the gap amplitude is several times larger than the size in the experimentally observed gap of the cuprate superconductors, if employ a widely accepted parameter mapping. The effective attractive interaction of carriers in this state is then estimated again to be extremely large in the order of or even larger than the nearest neighbor transfer energy t in the Hubbard model. We discuss implications of the results.
Finally in the supplementary information the reader can access to a fully functional open source code that was used to generate the data. The code is an extension based on the open source code mVMC [34].

A. Green's function
We present here a very general scheme to estimate onebody Green's function from the Lehman representation.
G e σ (k, ω) = Ω|ĉ kσ Hereĉ kσ (ĉ † kσ ) annihilates (creates) an electron of momentum k and spin σ. This approach requires the knowledge of the ground state |Ω of an HamiltonianĤ. The "hat" notation is used here to represent an operator as opposed to any matrix representation. Here η, is a small real number, a Lorentzian broadening factor.
In the ED, the Lehman representation can be evaluated explicitly since we have a complete and exact representation of the Hamiltonian eigenstates, but it is amenable only to small clusters. To evaluate this Green's function for the cases not amenable to the exact diagonalization, we can use a method similar to the approach used to calculate the spin and charge dynamical structure factors by exhausting important subspace of the Hilbert space for the excitations [24][25][26][27][28]. Namely, the idea of the present method is to restrict Hilbert space of the one particle or hole excited sector of the Hamiltonian to a set of vectors: n |Ω for electron excitations, (5) whereÂ n andB n are operators that together conserve the number of electrons N e and momentum k. Here |Ω is an approximate ground state of the N particle sector obtained by our variational Monte Carlo method for the ground state. Note that for the Krylov basis of excitation, used in ED,B n =Î andÂ n =Ĥ n where n is the number of band Lanczos iteration. Usually in ED, at every iterations of the band Lanczos method, the excited states basis is orthogonalized to every other excitations. But it is possible and sometime more convenient to work in the non-orthogonal basis.

B. Non-orthogonal basis
In the most general case, the excitated states chosen in our method given below are not orthogonal to one another. Then we introduce a number of overlap matrices: whereÎ is the identity operator and the matrix notation is expressed with the indices m and n. It is important to distinguish the operator notation (with a "hat" here) from the matrix notation since they are different in an non-orthogonal basis. Indeed, the matrices O k,mn are representations of the identity operator in this nonorthogonal basis.
Using the same basis, we evaluate the effective Hamiltonian matrices: In general, the basis set on the restricted Hilbert subspace is nonorthogonal and thus we need to solve the generalized eigenvalue problem within this subspace as: where H k and O k are matrices whose (m, n) components are given in Eqs. (7) and (6), respectively in the basis of |h kn or |e kn . The solution of this generalized eigenvalue problem is represented by the -th eigenvalue E k and its corresponding eigenstate coefficients of its eigenvector |E k represented in the basis of |h kn or |e kn . Note that the orthogonality of the eigenvector |E k is represented by E k |O k |E kj = δ ,j . The eigenvectors expand the basis of the subspace defined by the restricted Hilbert space determined by the choice of non-orthogonal excitations. We can insert the complete set of this subspace l |E kl E kl | in both Eqs. (2) and (3) to obtain: Note that |E e kl (|E h kl ) is the state in the N +1 (N −1) particle sector with momentum k. It is important to keep in mind that this is an approximation to the Green's function as we restricted the N +1 or N −1 particle sector by a variational form defined below in addition to the approximation to the ground state. By increasing the dimension of the excited subspace, we are able to systematically improve the representation of the excited states toward the exact one represented by the full Hilbert space. The accuracy can be tested by the convergence when we increase the number of variational states.
Note that we can adapt this formalism to re-express the band Lanczos algorithm. The main difference here is that we impose the translational symmetry from the beginning and we do not orthogonalize at each iteration aside from the variational form for the ground state. In fact, this is not an iterative technique, but for a chosen Krylov subspace, this formalism would give the same result as the band Lanczos method.
Our goal is to apply this formalism to obtain the spectral function of the strongly correlated fermion systems by combining with the variational Monte Carlo (VMC) method. An alternative derivation of Eqs. (9) and (10) can be found in Appendix A. This is the one implemented in the software provided in supplementary information.

C. VMC
In the variational Monte Carlo, we postulate an ansatz, a variational state |ψ that can be used to calculate the physical quantities associated with that state. To find a good approximation to the ground state, we optimize the variational state in order to minimize the energy measured [34,35].
The measurement of any operatorÂ can be done as: In this equation, the sum x |x x| is a complete set of every possible electronic configuration in the system. For large systems, however, it becomes computationally unfeasible to sum every configuration. Nevertheless we can still estimate this sum using a Monte Carlo sampling. The real-space configurations {x s } are generated with the probability ρ(x s ). Then, where N MC is the number of Monte Carlo sampling for the summation over x s , where s is the index to specify the particle configuration in the real space. For the ground state wavefunction, we employ the vari- Figure 1. Geometries of the system studied in this paper: a. is the one-dimensional lattice and b. is the two-dimensional square lattice. The excitation on site i is correlated with the presence of electron on neighboring sites. The charge of the excitation of typeĉ † iσn iσ |Ω is represented by the red arrow. The charges of the excitation of typeĉ † iσn i+δ,σni+δ ,σ |Ω is represented by the two blue arrows and the charges of the excitation of typeĉ † iσn i+δ,σni+δ ,σ |Ω is represented by the two green arrows.
ational state: where the variational parameter are f ij , g i and v ij . We define the variational ground state |Ω as the state |ψ that minimizes the variational energy E = ψ|H|ψ / ψ|ψ . To simultaneously optimize the variational parameters, the natural gradient method is applied. [36,37].

D. Dynamical VMC
Now, the excited states are constructed by applying the VMC calculation, because it is hard to generate excited states that reproduce the Krylov basis since the calculation ofĤ n is very expensive for n > 1 as every hopping of electron terms in the Hamiltonian produces a new Pfaffian evaluation. Instead we choose a basis wherê A n =Î andB n is a combination of different charge excitationsn iσ [28]. For example, the excitation basis we choose are |h iσn =ĉ iσ |ψ in and |e iσn =ĉ † iσ |ψ in where: |ψ in ={|Ω ,n iσ |Ω , n i+δn,σni+δ n ,σ |Ω ,n i+δn,σni+δ n ,σ |Ω }.
Here δ n and δ n are a combination of different neighbors to site i for each n. This choice is based on the physical reason, where the creation of an electron on site i is influenced by the presence of electron with both same and opposite spin on sites i + δ and i + δ respectively. This is illustrated in Fig. 1 for the two geometries studied in this article. Generally, we consider only excitation within a certain range δ = (δ x , δ y ) of neighborhood of the considered site i, noted as max(|δ x |, |δ y |) ≤ δ max where δ max is an integer that specifies the farthest neighbor considered in any direction. It is implicitly applied to both δ and δ . Under that threshold, we consider every combination of δ and δ that generate unique new excitation. The excitations have to be non-redundant (different). Otherwise this will lead to singular matrices for H k and O k . Note that |e iσn is in the electron-number sector with one-electron added to the ground state, while |h iσn is in the one-electron removed sector. Equation (18) is the simplest choice of basis to reasonably represent the essential part of low-energy excitation subspace of the Hamiltonian (as demonstrated in the result section). We calculate the average of ψ im |ĉ iσĉ † jσ |ψ jn and ψ im |ĉ iσĤĉ † jσ |ψ jn in the real space representation for every i, j combination (N 2 terms) after Monte Carlo sampling of the Markov chain and then Fourier transform (thanks to the translational symmetry): This produces a number of matrices as many as twice of the number of the momentum points matrices. We proceed with the same approach to evaluate the hole counterpart as well (O h k,mn and H h k,mn ) and solve the generalized eigenvalue problem (Eq. (8)) for both holes and electrons.
Together with Eqs. (1), (9) and (10), we estimate the Green's function from the VMC. We call the technique presented in this section dynamical VMC (dVMC). Appendix B discusses about details to optimize the computation speed of the calculation of these excitations.

A. Model
To examine the accuracy of the present dynamical VMC method, we show the benchmark test of the standard Hubbard model with the Hamiltonian on the 1D chain and 2D square lattice. The first term proportional to the transfer t is the kinetic energy and   the sum is restricted to the nearest neighbor pair. The second term represents the onsite Coulomb repulsion proportional to U . We impose the periodic boundary condition throughout this paper. In the presentation, the same color scale is used throughout this paper to make the comparison between different methods and size easier, except for Fig. 6. Unless specified, we use parameters U = 8, and the nearest neighbor transfer t = 1 as the energy unit. We add a broadening factor η = 0.2. The chemical potential (Fermi level) is determined so as to meet the occupied part of the integrated spectral weight with the number of electrons (N e = integer) given in our canonical ensemble simulation of VMC. Results obtained by the present method on lattices of N = 16 sites are first compared to ED results in order to benchmark the accuracy of the methods in one and two dimensions. The ED results were obtained using the open-source software HΦ [38] B. One-dimensional lattice In this section we test both the Mott insulator (N = N e ) and the doped Mott insulator in one-dimensional Hubbard model by comparing with the exact diagonalization results for the 16-site chain. The results at half filling are shown in Fig. 2(b) in comparison with exact diagonalization results in (a). We see that the present dVMC results show nearly perfect and quantitative agreement about the Mott insulating nature between the results of the present method and the ED in terms of the dispersion of the lower (occupied) Hubbard band below E F and the upper (unocccupied) Hubbard band above E F , in terms of their broadnesses, relative weights and the Mott gap sizes.
The results for the 64-site chain is shown in Fig. 2(c), which shows much higher momentum resolution. Both in 16-site and 64-site results, the lower (upper) Hubbard dispersion below (above) E F has split into two dispersions around the Γ (k = 0) and k = π points. The upper flat dispersion branch in the lower Hubbard dispersion was identified as the spinon branch and the lower steeper dispersion branch was identified as coming from the holon excitation in the analysis of the ARPES data for SrCuO 2 [39,40]. The present dVMC calculation correctly captures the spin charge separation in the 1D Hubbard model.
The doped case is shown in Fig. 3. The comparison between the dVMC result in (b) with the ED result in (a) for the doping concentration δ = 0.125 again shows nearly perfect agreement. The same doping for the 64site systems in Fig. 3(c) shows good agreement with the result obtained by Kohno using the dDMRG for 60-site system with the open boundary condition [41,42]. We still see both the holon and spinon bands. We also note the presence of the hole-pocket behavior appearing at the Fermi level [41]. We show the dependence on the number of excitations taken into account in Appendix C.

C. Two-dimensional square lattice
In this section, we examine the square lattice both for the Mott insulator (N = N e ) (see Fig. 4), and the doped Mott insulator (Fig. 5) again. Comparison of the dVMC results in Fig. 4(b) and Fig. 5(b) with the exact diagonalization results in Fig. 4(a) and Fig. 5(a), respectively, proves good accuracy of the dVMC even for the 2D lattice. The comparison of 8 × 8 lattice results in Fig. 4(c) at half filling and Fig. 5(c) at 12.5% hole doping with the result of the cluster perturbation theory [42,43] shows a fair overall agreement. Note that the results of the cluster perturbation has fewer momentum resolution and should not be taken as a sufficiently accurate reference. Therefore, a small discrepancy does not necessarily mean the insufficiency of the present result.
To have a better view of the superconducting physics, we reproduced the results on a bigger cluster. In Fig. 6, we examine a 12 × 12 cluster with the same parameters as Fig. 5(c). We see similarities between Fig. 6(a) and Fig. 5(c) but with much more details. To see the detailed structure of the gap, we first plot, in Fig. 6(b)-(d) the spectral weight in the Brillouin zone for energies at and near the Fermi level, namely ω = −0.3, 0.0 and 0.3 to estimate the position of the gap opening (in other words the locus of the minimum gap, which should form the Fermi surface when the gap closes). In principle, the locus of the gap opening can be inferred from the maximum intensity of A(k, ω = 0). The gap anisotropy can be seen along this trajectory. However, because of the limited discrete points in the Brillouin zone allowed for finite size studies, a better estimate of the gap opening position is estimated by an interpolation of A(k, ω) between neighboring momentum points using the data at small but nonzero ω as well. In fact, the rough estimate of the largest gap indicates that it appears at the antinode with the amplitude ∆ ∼ 0.3t. Therefore, the interpolation can be made efficiently by the choice of |ω| ≤ 0.3t. At ω/t = −0.3, 0.0 and 0.3, we also show thus obtained trajectory of the maximum A(k, ω) intensity by white dots and the corresponding interpolated spectral weight for these dots on the right, in panels (e), (f), (g). It is important to interpolate between different ω because the Fermi level determined from the consistency between the electron number in the occupied part of A(k, ω) and the given nominal number itself has uncertainty arising from the discrete k points and the broadening factor. In addition, because of the gap, it contains another uncertainty in identifying precisely the gap opening position in Fig. 6(c). Nevertheless, in Fig. 6(f), we can see clearly that the gap anisotropy is essentially expressed by the d-wave superconducting gap (∆(k) = ∆ 2 (cos k x − cos k y )), which has the maximum at the anti-nodes (k = (0, ±π), (±π, 0)) and closes at the nodes (k x = k y ), though the precise functional form of the gap is beyond the scope of the present paper.
Note that the number of excitation, for the case in Fig. 5(c) and Fig. 6, have been limited to 118 excitations, by keeping only the range 0 ≤ δ ( ) x ≤ 2 and 0 ≤ δ ( ) y ≤ 2 to suppress the statistical error within our allowed computational cost.
The ratio V d = ∆/F with F = P d (r → ∞) being the superconducting order parameter ∆ d is the measure of the effective attractive interaction to form the Cooper pair. From Fig. 9 in Appendix D, ∆ d is estimated to be 0.055. Then the attractive interaction is roughly estimated as V d ∼ 1.7t, which is extremely large, implying that the superconductivity in the Hubbard model is unrealistically strong, if we take the Hubbard model as a model of the cuprate superconductors with t ∼ 0.5 eV as employed in the literature. In fact, the gap size itself is estimated to be ∼ 0.3t, interpreted as ∼ 150 meV for the cuprates, which is several times larger than the gap amplitude around or less than 50 meV observed in the cuprates [1,2]. However, as we mentioned above, the true ground state is replaced by the charge inhomogeneous stripe state because of such a strong attraction.  ). The same is true for panels (c) and (f) and for panels (d) and (g), but in these cases we used linear interpolations to more precisely identify the paths away from the discrete k points. The superconducting d-wave gap structure is clearly seen on panel (f). Note that in this figure, (a) has the same color scale as the rest of the paper, but (b) to (g) has another color scale (shown in top right), in order to see more clearly the details in the superconducting gap and around the Fermi level.
A view for the origin of the effective attraction is the energy gain by the recovery of the electron coherence (fading out of Mottness), which grows in a nonlinear fashion with evolution of the doping. This nonlinear reduction of the kinetic energy forming a strong convex upward curve of the electronic energy as a function of the carrier density with the negative curvature signals the effective attractive interaction of careers, because the quadratic term obviously represents the effective electron-electron interaction [45]. The attractive interaction estimated from the negative curvature again has the same order consistently with the present value ∼ t [46]. By using the present method, we have shown that the d-wave gap of the Hubbard model has the energy scale of t (more precisely ∼ 0.3t), which is much larger than the value as the model of the cuprate superconductors. We note that more realistic understanding of the scale of the attraction of the cuprates and other available superconductors has to be reached by using the ab initio effective Hamiltonian that has realistic nonzero off-site Coulomb repulsions [44].

IV. SUMMARY AND OUTLOOK
In this paper we examined a newly proposed dVMC method to calculate the single-particle spectral function and the Green's function for strongly correlated electron systems. Although the proposed variational form of the excited states are simple and contains only one bare electron or hole added to the ground state dressed by composite operators diagonal in the particle number representation, the obtained spectral function rather accurately reproduces the exact structure in the benchmark test.
An application to the hole-doped square-lattice Hubbard model revealed that a d-wave superconducting state is induced by an effective carrier attraction, which is unexpectedly large in the order of the electron transfer resulting in much larger superconducting gap than that observed in the corresponding cuprates, if we study the charge uniform lowest-energy state as the ground state. It implies that the real cuprate superconductors have to be understood by taking account of more realistic factors such as the intersite Coulomb repulsion, which suppresses both the charge inhomogeneity and the superconducting order overestimated in the Hubbard model in comparison to the real existing superconductors.
Though it reached unprecedented and fruitful results, the obtained spectral function is not perfect with some discrepancy from the exact results. We note that spectral functions are very sensitive to the ground state |Ω used in the calculation, which are severely competing with other metastable states. The sensitivity requires a high accuracy of the ground state wave function, before calculating the dVMC Green's calculation. In addition, the form of the excited states has to be flexible enough particularly to represent low-energy excitations. Qualitatively different types of excitations ignored in the present work but presumably being important are the dressing by the spin-flip excitation such asĉ † i+δ,σĉ i+δ,σĉ † i,σ |Ω and by kinetic operators such asĉ † i+δ,σĉ j+δ,σ c † i,σ |Ω which are off-diagonal in the particle number representation. The inclusion of these excitations is an intriguing future issue to enhance the accuracy. Of course increasing the number of charge operators asn i,σni+δ n ,σni+δ n ,σ |Ω , n i,σni+δ n ,σni+δ n ,σ |Ω may also improve the accuracy.
In fact, the present dVMC method with such an improvement is expected to contribute to better understanding of the low-energy subtle structures such as pseudogap and effect of severe competitions among superconducting, charge and spin correlations of the doped Mott insulator.

V. ACKNOWLEDGMENTS
We acknowledge Youhei Yamaji, Kota Ido, Takahiro Ohgoe and Shiro Sakai for fruitful discussions. This work was supported by Fonds de recherche du Québec -Na-ture et technologies (FRQNT) and by a Grant-in-Aid for Scientific Research (No. 16H06345) from Ministry of Education, Culture, Sports, Science and Technology, Japan The authors are grateful to the MEXT HPCI Strategic Programs, and the Creation of New Functional Devices and High-Performance Materials to Support Next Generation Industries (CDMSI) for their financial support. We also acknowledge the support provided by the RIKEN Advanced Institute for Computational Science under the HPCI System Research project (Grants No. hp180170 and hp190145). A part of the computation was done at Supercomputer Center, Institute for Solid State Physics, University of Tokyo. and so on. These two optimizations greatly reduce noise in the resulting data.