Efficient world-line-based quantum Monte Carlo method without Hubbard-Stratonovich transformation.

By precisely writing down the matrix element of the local Boltzmann operator ([Formula: see text], where [Formula: see text] is the Hermitian conjugate pairs of off-diagonal operators), we have proposed a new path integral formulation for quantum field theory and developed a corresponding Monte Carlo algorithm. With the current formula, the Hubbard-Stratonovich transformation is not necessary, accordingly the determinant calculation is not needed, which can improve the computational efficiency. The results show that, the simulation time has the square-law scaling with system sizes, which is comparable with the usual first-principles calculations. The current formula also improves the accuracy of the Suzuki-Trotter decomposition. As an example, we have studied the one-dimensional half-filled Hubbard model at finite temperature. The obtained results are in excellent agreement with the known solutions. The new formula and Monte Carlo algorithm could be applied to various studies in future.

Introduction

Strongly correlated quantum many-body (SCQMB) systems possess a rich physical phenomenon, making them an essential topic in condensed matter physics. To overcome an exponentially increased dimensions of the Hilbert space, as well as the intrinsically strong correlation of SCQMB systems, scientists have strongly advocated for efficient and accurate numerical methods[1-7]. Over the last decades, several numerical methods have been proposed, such as the exact diagonalization (ED) method[8], the density matrix renormalization group (DMRG) method[9-11], as well as various quantum Monte Carlo (QMC) methods[12-14]. There are also many well-developed spin-based MC methods, which have been adopted for various correlated systems[15,16]. These numerical methods have played an important role in improving the understanding of SCQMB systems over the past few decades.

In ED methods, the Hamiltonian matrix of systems is directly diagonalized using advanced mathematical techniques, but they are limited to relatively small systems. The DMRG method shares some similarities with the ED methods but it can handle larger systems, making it a more effective method for low-dimensional systems. Some examples of commonly used QMC methods include determinant QMC (DQMC)[12,17,18], auxiliary field QMC (AFQMC)[13,19,20], and diagrammatic MC (DiagMC)[14,21]. Recently, we have developed a new method to directly calculate the ground state properties of elements of the Hamiltonian matrix, for a class of special systems[22,23].

The AFQMC method can be considered as a typical example of the DQMC method. For the DQMC method[12], a number of important advancements[24,25] and improved algorithms have been developed[17]. The DQMC method has been used to determine physical properties at zero[6,26,27] and finite temperatures[26,28]. The typical AFQMC method is a projection technique in which the operator continuously acts on a trial wave function; accordingly, the ground state properties[13,20,24,27] can be evaluated. To avoid the fermion sign problem or phase problem, AFQMC is usually implemented by an advanced restriction on the sample paths[29], and recently, the AFQMC has additionally been developed for finite temperatures[19,30-32]. DiagMC[14] combines the MC technique with Feynman diagrams of the perturbative expansion to calculate physical quantities[21,33-37], and several generations[38-41] and improvements[34,37,38,42,43] of this method have been developed based on the original DiagMC. For a more comprehensive introduction to various numerical methods, refer to review articles, such as those in Ref.[2,3].

In the DQMC methods, the quartic fermion operators (two-body term) are first decomposed into a quadratic one (one-body term) and a set of random external auxiliary fields by using the Hubbard–Stratonovich (HS) transformation[44-46]. The advantage of the HS transformation lies in that, in calculating the partition function, the trace can be easily evaluated by successive determinant calculations. However, the introduced auxiliary field, as well as determinants, increases the computation time heavily, which results in the computation time scaling cubic with the system size. In parallel to determinantal formulation mentioned above, the world-line formulation represents another catalog of QMC[47-51]. In the world-line QMC, the direct-space and imaginary-time is used to interpolate the representation of the fermion fields. The advantage of the world-line QMC is avoiding the time-consuming process of evaluating fermion determinants. However, in the world-line QMC, a closed path with non-zero weight is not always easy to sample in many-body wavefunction (WF) spaces. At present, several successful examples seem to be limited to a few specific Hamiltonians, and a universal algorithm is waiting to propose.

To further reduce the gap between the experimental studies and QMC calculations, more efficient numerical methods are needed to meet the current demands for studying SCQMB systems. Thus, there is an urgent need to either develop new numerical methods or optimize known ones. In this paper, we propose a new world-line QMC method by introducing a representation of the path integral formula in quantum field theory. This method can be used to calculate various properties of a system at finite temperature. The new formulation does not require the HS transformation, and not require determinants, therefore the calculation time can be significantly reduced. The results of the test calculation on the one-dimensional Hubbard model are in excellent agreement with exact values[52,53].

Proposed formula

To illustrate our method, we chose the simple but representative Hubbard model as an example. It is worth noting that our method can be directly extended to any model or Hamiltonian. The Hamiltonian of the Hubbard model reads:where denotes the creation (annihilation) of an electron with spin at the i-th lattice site. The first term on the right-hand side of Eq. (1) represents the one-body term, which accounts for the hopping of electrons between different sites, and is the hopping amplitude. The second term is the two-body on-site Coulomb interaction, where U represents the interaction strength, and is the number operator for spin at the site i. For convenience, the spin index () is omitted in the following description; yet it is included later on to prevent confusion.

One of the key steps of our method is to combine each off-diagonal term in the Hamiltonian and its Hermite conjugate into pairs, namely . Clearly is therefore a Hermitian operator. In the case of a general Hamiltonian, can be made up of the pair of a quartic fermion operator and its Hermitian conjugation. The purpose of this combination is to make as a Hermitian operator, and its eigenfunction can be easily obtained.

A many-body WF in the occupation number representation is labeled as . Here, the occupancy of the site i and j is explicitly given, while the occupancy of the rest of the sites is represented by K. For the site i and j, the occupation has four cases, which are labeled as . Here, indicates that there is no electron occupying the site i (j) site, while i (j) indicates an electron occupying the site i (j) site.

It is easy to prove that has only two eigenstates with non-zero eigenvalues:where the eigenvalues are equal to and , respectively. is the sign produced by the particle exchange as acts on . When an even number of exchanges occur, , otherwise, . Since and , we have . The remaining WFs orthogonal to and are also the eigenstates of , however, the corresponding eigenvalues are zero.

Since the operator and must share the common eigenvectors, the non-zero matrix elements of the local Boltzmann operator (LBO, , where can be any number) are only present in the following cases:where the remaining matrix elements of are equal to zero. Since the operator is diagonal, the non-zero matrix element reads:

Since , the matrix elements in Eqs. (3)–(4) are always positive, regardless of the value of . The matrix element in Eq. (2) is negative if is negative, otherwise, it is positive. Equation (2) produces the off-diagonal scattering in WF for the site i and j, but Eqs. (3)–(4) are the diagonal scattering in WF.

The partition function of a quantum many-body system is expressed as , where is the inverse temperature (or imaginary time) and refers the trace of in WFs space. According to the standard path integral formula, the imaginary time is divided into m time slices, where the partition function becomes , with the time step . Using the Suzuki–Trotter decomposition[44-46], the operator can be further decomposed as.

Finally, the partition function reads:

To calculate the partition function, a complete set of states are inserted between LBOs. For the purpose of our discussion, we number each LBO () in Eq. (6) from right to left as, 1st, 2nd, … i-th LBO. The WF following the k-th LBO is then denoted by the k-th WF. A closed world line (a closed WF sequence, or a closed path) in the WF space is labeled as , where is the s-th WF. is the associated Boltzmann weight of the world line .

Because only the matrix elements in Eqs. (2)–(5) are non-zero, for any closed world line, the Boltzmann weight has the following form:where correspond to the number of occurrences of the matrix elements in Eqs. (2), (3), and (5) for a closed world line , respectively. Evidently accounts for the number of occurrences of among , and determines the sign of the world line .

Equation (4) represents the weight of a world line based on the path integral formula in the current method. From this equation, one can find there are several novel features in our method: (1) The new formula does not include the HS transformation, thus it does not require the auxiliary field; (2) Our formula is not based on the determinant approach, thus the heavy calculation related to determinants is absent; (3) Our formula improves the accuracy of Suzuki–Trotter decomposition. According to the Suzuki–Trotter decomposition, the operator can be approximated by , where are any two operators. Considering the operator identity with . Thus, the error induced by the Suzuki–Trotter decomposition is directly related to . According to , the error in energies can be estimated by , where refers the ensemble average. If are Hermitian conjugate, then , which always contributes to the energy with an amount of . If are not Hermitian conjugate, as a first-order approximation, . Thus, we have enough reason to believe that, when are Hermitian conjugate, contributes the largest error in Suzuki–Trotter decomposition. Since in our method, the Hamiltonian is decomposed into Hermitian conjugate pairs, this error is automatically disappeared.

New Monte Carlo algorithm

To implement the new formula presented in Eqs. (6)–(7) into the QMC simulations, we have subsequently developed an efficient algorithm. Although it is easy to calculate the weighting of each path, it is relatively difficult to find a closed world line with a non-zero weight due to the fact that the weight of most paths is zero.

The choice of each WF is very significant for obtaining a closed world line with a non-zero weight. Here, we design an algorithm similar to the world-line algorithm[54] and the multiple time threading algorithm[55]. The current QMC algorithm contains two steps. To illustrate our method, we present an example in Fig. 1, in which the 4-site one-dimensional Hubbard model at the half-filled case is shown. Suppose is a closed world line from the last QMC step, the red line in Fig. 1 marks . The first step is to generate an intermediate world line () from a randomly selected WF in Specifically, a WF, say , is randomly selected from . In Fig. 1, the randomly chosen WF is marked by the arrow A. Then is scattered by the r-th LBO (), and a new WF is generated by . Next, the (r + 2)-th WF is generated by . This process continues until all the LBOs are cycled in the same sequence as in Eq. (6). In Fig. 1, this process starts from the arrow A to the right hand.

Figure 1

The graphic representation of the two-step sampling technique for a ring of four sites half-filled Hubbard model. The horizontal and vertical direction refers the time slices and lattice sites, respectively. “sij” is the shorthand of , indicating the scatting operator between two wavefunctions (vertical lines). The diagonal operator is omitted. In upper panel (step 1), the red line denotes a closed world line (), while the blue line denotes an open world line (). The arrow A and B mark intersections between and . In lower panel (step 2), a closed world line () is constructed by replacing the segment between arrows A and B in with that in .

In the above scattering process, does not produce any bifurcation because it is diagonal. However, for the operator , if one of the i-th or j-th sites is occupied and the other is empty, the scattering will be bifurcated. One side of the bifurcation corresponds to Eq. (2), in which the occupancy of the i-th and j-th sites is exchanged before and after scattering (corresponding to the diagonal line in Fig. 1). While the other path corresponds to Eq. (3), in which the wave function is unchanged before and after scattering (corresponding to the horizontal line in Fig. 1). For the bifurcation, we use a similar heat-bath algorithm[56] to select the new path. Assuming the first path following Eq. (2) with the probability of , and the second one following Eq. (3) with the probability of . Evidently . Note is different from that in Ref.[54].

After the scattering process finished, the intermediate world line is successfully generated. The blue line in Fig. 1 is the intermediate world line generated by the above scattering process. It needs to point out that, may not be a closed world line. In fact, in most cases, it is an open world line. In Fig. 1, is opened at the arrow A. The key point is that may have multiple intersections with . The intersection means at which the WF is identical in both and . In Fig. 1, the arrow A and B mark the two intersections. In the second step, the fragment between two randomly chosen intersections in is used to replace the corresponding part in , then a new closed world line () is constructed, illustrating in the lower panel of Fig. 1.

The acceptance rate of the new world line is determined by the ratio of two factors, namely , where is the Boltzmann weight of the new and old world lines according to Eq. (7). And is the total weight attached to the heat-bath sampling equal to . The product contains all the contribution of each bifurcation in the fragment of used in the new closed world line . should be deducted from the acceptance rate. If the change is accepted, the updated WFs are implemented. Otherwise, the unchanged WFs are implemented.

There are a few differences between the current method and previous ones[48,54]. (1) Except for the initially selected WF from , the scattering process is irrelevant to the rest WFs in , which is different from Ref.[54]. (2) In the current method, the sequence of in Eq. (6) can be arranged in any way, the only requirement is to combine each off-diagonal term in the Hamiltonian and its Hermite conjugate into pairs; (3) In the current method, there are two steps to generate a new closed world line (). The first step is to generate an intermediate world line () from a exist closed world line () by scattering process. The second step is to construct from and . The current procedure does not care whether is closed or not, but has at least two intersections with . In previous method[54], the scattering process is aimed to directly generate a closed world line in a single step. Because of this requirement, the previous method[54] usually needs a specific break-up or rearrangement of Hamiltonian. In comparison, the current algorithm is more straight forward than previous ones, and can be easily extended to any Hamiltonian. (4) As first feeling, one may expect there should be few intersections between and . However, the probability of finding a new closed world line using the current method is remarkably high, i.e., close to 100%. This may derive from the contribution of the Hermite pairs used in our method. (5) Because the scattering keeps the number of particles unchanged, similar to Ref.[54], the current method also works in a canonical ensemble, which is different from Ref.[48].

Tests on Hubbard model

The Hubbard model is a representative for studying typical SCQMB systems. By varying the strength of U, the Hubbard model can describe different systems from the weakly coupled case to strongly coupled case. Usually U = 8, 4 and 2 correspond to the strongly coupled, the intermediate coupled, and the weak coupled cases, respectively[57]. For both the strongly and intermediate coupled cases, the mean field method fails. In this paper, our major calculations will be concentrated on U = 8 and 4.

Many studies based on the Hubbard model have been carried out to investigate the metal–insulator transition, superconductivity, and magnetic properties caused by electronic correlation. Lieb and Wu[58] have obtained the exact solution of the ground state for the half-filled one-dimensional Hubbard model. In the last few decades, various theoretical calculations have been carried out to study the one-dimensional Hubbard model. However, few studies have been conducted using this model at finite temperatures[52,53,59-64].

To illustrate the reliability of our new formula and the new QMC algorithm, we have studied the one-dimensional half-filled Hubbard model at finite temperature. In all calculations, is taken to have units of energy and is set to = 1.0. The one-dimensional Hubbard model with the number of lattice size of N = 6, 12, and 24 have been systematically studied. For each system, the strength of the interaction has also been investigated for U = 2, 4, 6 and 8. To determine how the simulation time scaling with N, we also calculate a few larger systems with N up to 96. For most simulations, the total number of QMC steps at each temperature or m is more than 107, where the first third of the steps are used to equilibrate the system and the remaining two thirds of the steps are used to calculate the physical properties.

It needs to point out that, like most QMC methods, our new formula could not give a general solution for the sign problem too. However, for particular special cases, the sign problem is not encountered[26,65], for example, in one-dimensional Hubbard model. By choosing appropriate boundary conditions (periodic or antiperiodic) in one-dimensional Hubbard model, in Eq. (2) can be always positive, thus the sign problem can be avoided in current method. This is why we choose the one-dimensional Hubbard model.

The energy, double occupancy, local magnetic moment, and spin correlation functions have been calculated in the temperature range of 0.05 to 4.0. According to thermodynamics, the energy of systems can be calculated as . The double occupancy, denoting the probability of two electrons occupying one site, is written as . is the local magnetic moment, where denotes the z-component spin operator at the i-th site. The nearest-neighbor and next-nearest-neighbor spin correlation functions are defined as and .

The convergence test on is shown in Fig. 2 for N = 6, where the upper and lower panels present the data for T = 0.25 and 0.5, respectively. It can be seen that, with an increase of the number of time slices (m), the energy converges quickly. For T = 0.25, as m = 80, corresponding to τ = 0.05, the QMC results approach the exact value with a difference of approximately 1% and 5% for U = 4 and 8, respectively. Owing to the intrinsic characteristics of the path integral formula, the exact value can be obtained only as τ approaching to zero. To estimate the QMC energy at τ = 0, the extrapolation to τ = 0 is performed. For this purpose, the data shown in Fig. 2 is fitted by , where are fitting parameters. We find the value of is around 1.5, which is weakly dependent on U and slightly reduces with the increase of N. is the extrapolated energy at . The exact energy () and extrapolated energy () are shown in Fig. 2. The results show that our QMC method does converge to the exact value at .

Figure 2

Energy per site as a function of the number of time slices (m) for systems with six lattice sites at temperature of 0.25 (upper panel) and 0.5 (lower panel) at both U = 4 and 8. The symbol, dashed lines, and solid line represent the QMC data, the exact value, and the fitting curves, respectively. and refer the QMC energy extrapolated at and the exact energy, respectively.

We have tested the convergence speed for several systems with different N and U. To compare the convergence speed, we define a convergence criterion (), at which the derivative of energy with m is less than 0.0005. This criterion is equivalent to that the change in energy is less than 0.0005 as m increases by a unit. Figure 3 shows how changing with U and N. It can be seen that the convergence speed does not change significantly with N, but decreases with the increase of U. For most cases, τ = 0.05 is already a good approximation. In the following QMC simulations, τ is fixed at the value of 0.05.

Figure 3

The convergence criterion () varying with the on-site Coulomb interaction strength (U). Square, triangle and circle symbols represent the systems with lattice size of N = 6, 12 and 24, respectively.

Since our QMC method is independent of determinant, the simulation time should have much better scaling with N and m. To check this point, we have calculated the simulation time as a function of N and m at fixed U = 4. These calculations are performed on a desktop computer with the CPU basic frequency of 3.20 GHz (Intel Core I7-8700) and a serial QMC program. The simulation time for 200 thousand MC steps is calculated for various systems, which is summarized in Fig. 4. From this figure, one can see that, the simulation time has the linear scaling with m (upper panel of Fig. 4) and the square-law scaling with N (lower panel of Fig. 4). This computational cost is far lower than other QMC methods involving determinant calculations, and is comparable with the common first-principles calculations. It should be stressed that, our current QMC code can be further improved for the higher efficiency.

Figure 4

The simulation time as a function of the number of time slices m for N = 6 (upper panel) and the number of lattice size for m = 40 (lower panel). The symbol and solid line represent the QMC data and the fitting curves, respectively. Here the linear fitting to QMC data is adopted in upper panel, while the quadratic function is used to fit QMC data in lower panel.

Although the Hubbard model is a benchmark system for testing various QMC methods, a single model may be not enough to demonstrate the advantage of our method. To remedy this issue, we have done an analysis of the computational complexity for our method. The amount of calculation mainly consists of two parts: (1) finding a new closed path. This part needs to calculate the scattering of all LBOs to adjoint WFs, which needs Nm operations; And after each scattering, the comparison between the new and old WFs is preformed, which needs N operations. Thus, the total amount of calculations is scaling as Nm. (2) The calculation of Eq. (7). In this step, the scattering matrix of each LBO is calculated, the corresponding computational costs is proportional to Nm too. Combining these two parts, the total amount of calculation scales as N2m. This is also consistent with our test results (Fig. 4).

In the following, we will present detailed calculations of various physics quantities for the system with six lattice sites. Figure 5 depicts the energy, double occupancy, and specific heat via temperature for U = 4 (left panel) and 8 (right panel) of the system with six lattice sites. The energy calculated by the QMC simulation is in excellent agreement with the exact value for the entire range of temperatures within the error bar (upper panel of Fig. 5). Compared to U = 4, there is an evident plateau in the temperature range of 1.0 to 2.0 for U = 8. The plateau reflects the fact that the on-site Coulomb interaction has a strong effect in suppressing the occurrence of double occupancy. The change in double occupancy with temperature in middle panel of Fig. 5 supports this conclusion.

Figure 5

Energy (upper panels), double occupancy (middle panels), and specific heat (bottom panels) as a function of temperature for U = 4 (left panels) and U = 8 (right panels) of the system with six lattice sites. The QMC results (symbols) are in excellent agreement with the exact value (dashed lines) across the entire range of temperatures.

From the middle panel of Fig. 5, it can be seen that, from the high temperature to the lower temperature, the double occupancy first decreases and then increases for U = 4 and U = 8. Although the double occupancy increases at low temperatures, the total energy is further reduced with a corresponding decrease in temperature. This reflects the fact that a small increase of the delocalization doublon further decreases the total energy[52,66]. The minimum indicates the degree of localization of electrons is the largest, which corresponds to the maximum in local magnetic moment in upper panel of Fig. 6. For U = 4, the increase in the double occupancy is more evident than when U = 8 at lower temperatures. At a fixed temperature, the double occupancy is larger than it is for U = 8, demonstrating how the on-site interaction has a noticeable impact on the formation of the double occupations.

Figure 6

Local moment () and spin correlation functions () as a function of temperature for U = 4 (left panel) and U = 8 (right panel) of the system with six lattice sites. , , and represent the local magnetic moment, the nearest-neighbor spin correlation, and the next-nearest-neighbor spin correlation, respectively.

The specific heat as a function of temperature is shown in the lower panel of Fig. 5. To calculate the specific heat, the exponential fitting method[67-69] is adopted with the fitting form of , where , and are the fitting parameters. In this study, the value of M was 8. From Fig. 5, it can be seen that there are two obvious peaks in the specific heat for U = 8. Specifically, there is a narrow peak at low temperatures and a broad peak at high temperatures. In contrast, for U = 4 the two peaks become much closer and begin to merge together. The stronger interaction strength U, the more obvious the peak. The structure of the obtained specific heat peak is consistent with previous findings[52,61,62]. It is believed that this feature in the specific heat is associated with the spin-wave excitations at low temperatures and the single-particle excitations at high temperatures. Thus, spin fluctuations and charge fluctuations are dominant at low and high temperatures, respectively, which can be highlighted by the correlation functions. The trends of spin correlations in low panel in Fig. 6 have been correlated with the peaks in specific heat. The results we have obtained are consistent with known values.

Figure 6 shows the local moment ( and spin correlation functions ( as a function of temperature for U = 4 (left panel) and U = 8 (right panel) of the system with six lattice sites. With an increasing temperature, reaches its maximum at a certain temperature and then gradually decreases. The maximum value obtained for indicates that at this temperature, the degree of localization of electrons is the largest. The degree of delocalization of electrons reflects the formation of doublons; therefore, the trends of local moment and double occupancy are reversed, as shown in the upper panels of Fig. 6 and the middle panels of Fig. 5. From the lower panel of Fig. 6, one can see that is less than zero, which indicates an antiferromagnetic order at a finite temperature, which leads to the emergence of a specific heat peak at a lower temperature. As the temperature increases, decreases and tends to zero, reflecting a weakened antiferromagnetic order; this is in agreement with the exact results. In contrast, is greater than zero and gradually reduces to zero with a corresponding increase in temperature. It can be seen that, the error bar in for U = 8 is relatively large, which may be due to the limited simulation time, or the intrinsic large fluctuations in spin correlation functions. Fortunately, the trends of are in general consistent with the results of Shiba[52,53].

Summary

We have proposed a path integral formula in field theory and a corresponding world-line quantum Monte Carlo algorithm. The remarkable feature of the current method is that neither determinants nor the HS transformation is needed, which does strongly improve the accuracy and efficiency of Monte Carlo simulations. As an example, we have calculated the thermodynamic quantities and correlation functions of the one-dimensional Hubbard model at finite temperature. Our results are in excellent agreement with the exact values, confirming the reliability of our method. The most encouraging thing is that the computational cost has the square-law scaling with the size of systems. We believe that the current approach could be widely used in future.