Importance-sampling FCIQMC: Solving weak sign-problem systems

The fermionic sign problem is one of the most fundamental limitations of quantum Monte Carlo (QMC) methods.^1,2 In projector QMC methods, where the aim is to project out a specific state—usually the ground state—through a stochastic evolution of signed walkers, the problem manifests itself through the exponential growth of an undesired distribution of walkers of “bosonized” or “stoquastized” character, leaving the desired fermionic distribution lost in the noise of the simulation. This instability leads to diverging variances in expectation values of measured observables.

There are three known strategies to counter this problem. The first is to place constraints on the evolving distribution of walkers to prevent the growth of the undesired distribution. The fixed-node approximation of Diffusion Monte Carlo (DMC),^3–5 the constrained path and phaseless approximations of auxiliary-field quantum Monte Carlo (AFQMC),^6–12 and the initiator approximation¹³ of the full configuration interaction quantum Monte Carlo (FCIQMC) method¹⁴ all fall into this category. The difficulty is that these approximations generally lead to biased results (unless an exact form of the constraint, e.g., via a symmetry, can be applied), and furthermore, these approximations are either uncontrolled or can only be controlled (reduced) through time-consuming protocols, such as multi-determinant trial wavefunctions in the fixed-node and phaseless approximations, or by increasing the walker number in the initiator approximation. These methods are ultimately scaling exponentially in the computational resources needed to achieve a desired accuracy. The second strategy is to perform unconstrained projections, attempting to measure a transient fermionic signal before the undesired solution becomes dominant. Such simulations are practical only for very small systems and otherwise cost exponential time. The third strategy, which is advocated in the FCIQMC method (without the initiator approximation), is to introduce walker annihilation processes into the algorithm, which destabilize the undesired solution¹⁵ and allow the desired fermionic solution to emerge without bias. The difficulty here is that a minimum number of walkers is required to achieve this, and this minimum number is generally found to scale with the size of the Hilbert space. This makes the algorithm exponentially costly in memory, though usually with a small prefactor. When exact FCIQMC is applicable, the quality of the sampled wavefunction is comparable to that obtained by exact diagonalization, but with a lower memory requirement.

When FCIQMC is used with a small number of walkers, the unphysical stoquastic solution becomes dominant and the fermionic solution is lost in noise. The stability of the fermionic solution can be quantified by the energy gap between the exact fermionic ground-state and that of the undesired solution.¹⁵ However, as we will show in this paper, the stoquastized gap alone does not fully describe the computational difficulty to overcome the sign problem. In particular, we will show that two Hamiltonians with exactly the same gap can, nevertheless, require very different minimal numbers of walkers to converge. We show that guiding wavefunctions can be used to further mitigate the sign problem by reducing the minimum of walkers required to fully resolve the sign problem, thereby arriving at an exact scheme with lower computational cost. (Previously, it has been found that by employing a trial wavefunction, a “fixed-node Hamiltonian” can be defined by which this stoquastized gap can be reduced at the expense of introducing an approximation to the Hamiltonian.)^16–19 

If a Hamiltonian in a particular basis has no sign problem, then its stoquastized gap is identically zero; conversely, if this gap is found to be zero, then that Hamiltonian has no sign problem. A weak sign problem is one in which the gap is non-zero, but small. There are some systems, such as periodic 1D Hubbard models, in which the gap remains small even as the system size grows (i.e., it is not size-extensive). In these cases, the FCIQMC algorithm can be applied to very large systems, essentially accessing the thermodynamic limit. In other systems, the gap grows with system size, ultimately limiting the system sizes, which can be treated. Even here, one can still obtain exact-diagonalization quality results for systems far exceeding conventional exact diagonalization.

A characteristic of weak sign-problem fermionic systems is that their wavefunctions tend to be highly spread out throughout the Hilbert space, somewhat similar to the related stoquastized systems (in which all imaginary-time Feynman path contributions are strictly positive, leading to diffuse, non-compact wavefunctions). For such systems, methods based on truncations of the Hilbert space will lead to large errors. Therefore, techniques such as the initiator approximation of FCIQMC, or similar selected CI methods of quantum chemistry, struggle to describe such systems. Paradoxically, therefore, weak sign-problem systems are highly challenging systems, and it is necessary to devise alternative techniques—not based on truncations—to address them. The FCIQMC method without the initiator approximation is one such method, but even there, it needs to be adapted to efficiently yet faithfully sample the exact fermionic ground-state. Devising such a method is the main purpose of the present study.

It is an interesting fact that not every fermionic quantum system has a sign problem in FCIQMC, even in non-trivial cases. In some systems, it turns out that it is possible to subdivide the Hilbert space of the problem in such a way that, during the propagation step of the FCIQMC algorithm, all walkers that arrive at a given point in the Hilbert space always have the same sign, even though there are both types of walkers distributed throughout the space. Such a case occurs in one-dimensional Hubbard models with open boundary conditions, irrespective of the filling, as long as one uses a real-space basis.²⁰ 

The Hamiltonian for a Hubbard chain of length ℓ is given by

$ĉi,σ†$ and $ĉi,σ$ are the fermionic creation and annihilation operators at lattice site i with electron spin σ, respectively. $n̂i,σ=ĉi,σ†ĉi,σ$ are the corresponding number operators. t is the hopping amplitude, and U is the on-site interaction. The Hubbard model is a seemingly simple lattice model but yet already shows the rich physical properties of more realistic system setups^21–23 and has been studied with a variety of different computational methods.²⁴ 

For open boundary conditions—sites ℓ and 1 are not regarded as neighbors—the system is sign-problem-free for all fillings. For periodic (anti-periodic) boundary conditions—sites ℓ and 1 are connected with hopping amplitude −t (t)—there is no sign problem for an odd (even) number of up-electrons and an odd (even) number of down-electrons.

In these systems, calculations of systems with up to 150 sites in the case of the one-dimensional Hubbard chains have been shown to be calculable with FCIQMC.²⁰ In such large spaces, another systematic bias called population control bias, which is independent of the sign problem, becomes notable. It has been shown that the application of importance sampling in combination with an a posteriori reweighing procedure can practically remove the bias.

When moving away from the filling rules in the periodic 1D Hubbard model mentioned previously or moving to 2D Hubbard systems, there is an inevitable sign problem. However, we will show that these systems still show only a weak sign problem when using a real-space basis, meaning that the number of FCIQMC walkers to fully resolve the sign problem is very small compared to the total Hilbert space size.

To illustrate that in a real-space basis the sign problem of the Hubbard model is comparatively weak, we calculated ΔE^stoq of two half-filled Hubbard 4 × 4 systems. For U/t = 8, in a real-space basis, we obtain ΔE^stoq = 3.792. In a reciprocal-space basis, we get 97.201, which means an increase of almost two orders of magnitude. These numbers are also shown in Fig. 2. Even for U/t = 4—where the k-space representation typically is much more suitable²⁵—ΔE^stoq in the real-space basis is still much lower than in the reciprocal-space basis (8.401 compared to 40.720).

The sign problem is especially weak if the lattice geometry still resembles the geometry of a sign-problem-free 1D system, such as a two-legged ladder. Still, the minimum number of walkers may be too high to be affordable.

In this paper, we analyze the FCIQMC algorithm applied to weak sign-problem fermionic systems, epitomized by Hubbard models in a real-space basis in periodic 1D and ladder geometry (see Fig. 1). In such problems, the initiator approximation is severe, while the exact algorithm can be applied affordably using importance sampling. We observe a significant reduction of the minimum number of walkers in weak-sign-problem systems. We use the method to report unbiased results for the interaction regime of U/t = 8 at and close to half-filling.

FCIQMC is an imaginary-time projection method to (in its simplest formulation) sample the ground state wavefunction |Ψ₀⟩ of a Hamiltonian $Ĥ$ by solving the stationary Schrödinger equation.¹⁴ The master equation is given by

with S being a scalar parameter, the shift. Starting from an initial wavefunction |Ψ(0)⟩ that has non-zero overlap with the true ground state, |Ψ(τ)⟩ will converge to the ground state in the large-τ limit.

FCIQMC applies this master equation in a finite basis of Slater determinants (although recently FCIQMC has been extended to a spin-adapted version working in a basis of configuration state functions^26,27). The evolving wavefunction coefficients are represented by a number of walkers N_i that reside on a specific Slater determinant |D_i⟩. $Ĥ$ is, then, stochastically applied to the instantaneous walker distribution in discretized time steps Δτ as

The stochastic application of $Ĥ$ consists of three steps: spawning, death/cloning, and annihilation. In the spawning step, each walker attempts to spawn from the determinant |D_i⟩ it is sitting on to a determinant |D_j⟩ that it is connected to by a non-zero Hamiltonian matrix element H_ji. The spawning attempt succeeds with a probability $pspawn=ΔτHji/pjigen$⁠, where $pjigen$ is the probability that |D_j⟩ is selected as a target determinant starting from |D_i⟩. In the death/cloning step, the walker population N_i of every occupied determinant is multiplied by Δτ(H_ii − S). If after this step the walker population is lower than a certain minimum occupation threshold (mostly chosen to be 1), the population is stochastically rounded either up to the occupation threshold or down to 0. If it is rounded to 0, the determinant is removed from memory until a new spawn happens onto this particular determinant. This makes the instantaneous representation of the FCIQMC wavefunction very memory-efficient. The annihilation step is an important step for reducing the sign problem. As the last step of an iteration, all spawns to the same target determinant are communicated to the same processor. There, positive and negative spawns cancel each other. This is easily possible since FCIQMC works in a discrete basis of Slater determinants unlike, e.g., diffusion Monte Carlo (DMC) where such immediate cancellations are not straightforwardly possible and, therefore, of limited effectiveness.⁵ 

The ground-state energy can be measured by averaging either the projected energy or the shift. The projected energy is defined as

where |D₀⟩ is the reference determinant (which in real-space Hubbard systems close to half-filling is always a Néel determinant). The shift is not only used to control the population but can also be used as an energy estimator. It is updated according to

every A iterations. The third term is an extension to the usual shift update formula and was presented in Ref. 28. The first logarithmic term just leads to a constant population. The second term has the desired population N_w in the denominator so that the average total population will converge to that value precisely. The prefactor γ²/4 ensures critical damping in a simplified model and still proves useful in actual simulations.

Typically, in ab initio systems, FCIQMC is applied in conjunction with the initiator approximation¹³ and has been successfully applied in various ab initio systems,^29–32 as well as the Hubbard model in a reciprocal-space basis for U/t ≤ 4.^25,33 In the initiator method, determinants are dynamically marked as initiators during the simulation when their walker population exceeds the preset initiator threshold t_init. Only initiators are allowed to spawn to every other determinant since they are deemed to be sign-coherent because of their well-established population of a single sign. Non-initiators can only passively acquire population and spawn to occupied determinants only, with the rare exception that two simultaneous non-initiator spawns on an empty determinant are also allowed. This leads to an effective, dynamically updated truncation of the Hamiltonian. This introduces a bias into the sampled wavefunction and energy estimators. The initiator approximation is not suitable for treating the real-space Hubbard model for U/t ≲ 16 where the wavefunction is strongly spread-out across the Hilbert space. Even for t_init close or equal to the occupation threshold, the initiator truncation becomes so severe that the energy estimates are very slowly converging with respect to N_w, even when using improved formulation of the initiator method, such as the adaptive-shift method.^34,35 Therefore, in the following, we will discuss calculations and algorithmic improvement within non-initiator FCIQMC only.

In the absence of any annihilations, the FCIQMC wavefunction would collapse to the ground state of the so-called “stoquastized” version $Ĥstoq$ of the original Hamiltonian.¹⁵ The matrix elements of the stoquastized Hamiltonian are defined as

In other words, the matrix elements of $Ĥstoq$ have the same magnitude as $Ĥ$⁠, but the matrix is a sign-free version of it. It is easy to see that the lowest eigenvalue $E0stoq$ of $Ĥstoq$ is always lower than or equal to the lowest eigenvalue E₀ of $Ĥ$⁠.³⁶ The larger the stoquastized gap $ΔEstoq=E0−E0stoq$⁠, the larger is the tendency of the shift to collapse to the ground state with eigenvalue $E0stoq$ instead of the sought-for true fermionic ground state with eigenvalue E₀. On the same note, the overlap with any reasonable trial wavefunction vanishes because the fermionic solution is exponentially suppressed. This also makes the determination of the ground-state energy using projected energy estimators impossible.

ΔE^stoq is an important concept to quantify the sign problem. It can be regarded as a global quantity that quantifies how far the stoquastized and fermionic differ. However, there is no monotonic relation between the stoquastized gap and the actual computational effort (i.e., the minimum number of walkers required) to solve a system with FCIQMC. In a practical simulation, the actual sampling of the wavefunction is also important. This fact will be demonstrated and exploited later in this paper.

To render the correct fermionic solution of a given sign-problematic Hamiltonian with FCIQMC, the Hilbert space has to be occupied with a minimum number of walkers simultaneously. This is such that the necessary annihilations can actually take place and the fermionic solution can emerge. Still, the minimum number of walkers N_min to get the numerically exact ground-state solution and overcome the sign problem is typically smaller than the number of determinants $|H|$ in the Hilbert space. The minimum number of walkers is also called the annihilation plateau since in most cases, it can be identified by a plateau in the walker number in the walker growth phase, caused by walker annihilations. However, especially in the real-space Hubbard systems studied in this paper, the plateau can be hard to identify. Therefore, we use a different method to determine N_min, which we call the fixed-N₀ method. In this method, the population on the reference determinant is held fixed, rather than the total population. In order to achieve this, the shift is updated according to

N₀(τ) is the population of the reference determinant at imaginary time τ and j runs over all connected determinants, including j = 0. The shift is now equivalent to E_proj, and there is only a single energy estimator. Plugging Eq. (7) into Eq. (3) leads to

so the population on the reference determinant is held exactly constant at a specified value N₀ and especially at a given sign. This ensures sign coherence automatically at the expense of a total walker population N_w that is not predefined. For N₀ = 1, the averaged total population $Nw̄$ would converge to N_min exactly. Choosing N₀, this low, typically, leads to large fluctuations of N_w, however, making it difficult to average. Therefore, in all calculations, in this paper, we will use N₀ = 50, which will lead to $Nw̄$ slightly higher than N_min.³⁷ 

As mentioned in the Introduction, the 1D Hubbard model is sign-problem-free for open boundary conditions for all fillings and for periodic boundary conditions for an odd number of up-electrons and odd number of down-electrons. A proof for this is given in  Appendix A.

Moving away from these system setups to different fillings or going to 2D will inevitably introduce a sign problem. The strength of the sign problem, however, strongly depends on the basis representation, the geometry of the lattice, and the filling.

As can be seen in Fig. 2, the sign problem in a real-space basis is typically much weaker than comparable systems in a reciprocal space basis, i.e., ΔE^stoq is an order of magnitude smaller.

We next analyze the sign problem of 1D systems with even n_↑ or n_↓ with periodic boundary conditions. According to Eqs. (A3) and (A4), there are non-negative matrix elements in the Hamiltonian; however, the sign problem is non-size-extensive. This means that the sign problem is actually weakening with increasing system size.

To understands this, let us see how the sign problem emerges. Let $Pji=⟨Dj|1−ΔτĤ|Di⟩$ be the amplitude to propagate a walker on |D_i⟩ to |D_j⟩ in one time-step of Δτ. Then, the amplitude to propagate a walker from |D_i⟩ to |D_j⟩ after K such steps is given by the sum over all possible Feynman paths of length K, which start at i and end at j,

For ground-state problems, we are interested in this quantity in the K → ∞ limit, and we will refer to this transition amplitude as

Given the stoquastized problem, the amplitude to propagate a walker from |D_i⟩ to |D_j⟩ is given by

with the corresponding ground-state quantity being

The stoquastized quantities $Pji(K),stoq$ and $Cjistoq$ are by definition non-negative. A sign problem, i.e., a gap between the fermionic and stoquastized ground states, is created when there are at least two paths in $Pji(K)$ with different signs, leading to

The maximal sign problem is created when there is full destructive interference in the fermionic contribution, i.e., C_ji = 0. Let us now consider what happens in different geometries.

In the 1D case, the number of pathways between two determinants is severely limited due to the lattice geometry. In open boundary conditions, there is only one essential pathway—thus no sign problem—in periodic boundary conditions, there are only two essential pathways that could destructively interfere. The sign problem measure $Cjistoq−Cji≥0$ is reduced if the chain length ℓ is increased or if the on-site interaction U/t is increased. An illustrative explanation for these relations is given in  Appendix B. This leads to a non-extensive sign problem for 1D Hubbard systems, which is shown in Fig. 3(a) in terms of the stoquastic gap as a function of the chain length. For systems with over 100 sites, the sign problem becomes negligible even off-half-filling and with periodic boundary conditions as $E0stoq$⁠, then, becomes equal to E₀ within statistical errors. Thus, systems of this size can be solved with less than 10⁸ walkers when accounting for the population control bias.

On the same note, this leads to the counterintuitive situation that for sign-problematic 1D systems of intermediate size (e.g., around 50 sites), where $E0stoq$ is still smaller than E₀ outside of statistical errors, is much harder or even impossible to solve. Even if the non-size-extensive gap is already very small at this system size, due to the very large Hilbert space size, a significant or even unreachable number of walkers would be needed to establish a sign-coherent solution.

In the 2D case for a square lattice, the length of the elementary sign-problematic cycles does not change with the increasing system size. It always consists of a four-site plaquette. Therefore, the special circumstances that cause the non-extensivity of the sign problem in the 1D case no longer exist here. In addition to the elementary cycles, adding more four-site plaquettes even increases the number of sign-problematic cycles. This leads to an extensive sign problem, which means that the stoquastic gaps increase with system size. This is shown in Fig. 3(b). This severely limits the treatable system size. Depending on the exact shape of the lattice, however, the sign problem can—even though it is size-extensive—still be weak and systems well beyond the capability of exact diagonalization can be treated numerically exactly when FCIQMC is enhanced by an importance sampling scheme using a simple guiding wavefunction. An example for a Hubbard system class with a very weak sign problem is the ladder systems, shown in Fig. 1, which are 2 × ℓ Hubbard lattices with a 2 × 2 primitive cell. Ladder systems have been studied mostly using the density-matrix renormalization group (DMRG) as they still show a relatively low entanglement entropy, making it possible to calculate larger clusters.³⁸ 

We end this section by noting that our finding that the real-space representation results in a weaker sign problem compared to the momentum space representation is in line with a similar conclusion arrived at by Blunt, in a study of fixed-node and partial-node FCIQMC,¹⁹ for ab initio systems. There, he found that fully localized basis sets tend to lead to weaker sign problems (as measured by the population plateau), compared to split-localized and canonical Hartree–Fock or natural orbitals. Taken together with our results, we believe this finding may be quite general, namely, real-space or localized basis sets—while leading to highly spread-out, non-compact, wavefunctions—also lead to the weakest sign problems. Put another way, we hypothesize that orbital representations that lead to compact wavefunctions tend also to lead to severe sign problems, implying that the observed compactness may be due to large amounts of annihilation.

We support this by the following argument: Suppose that we perform a deterministic simulation using the FCIQMC master equation [Eq. (3)]. At equilibrium, the net contribution onto |D_i⟩ is ΔN_i(τ) = 0. When all other N_j equal their respective true C_j coefficient, N_i is then given by

N_i will be small if (i) H_ii is large, (ii) all connected N_j or their respective matrix elements H_ji are small in magnitude, or (iii) there are a lot of canceling contributions from connected determinants |D_j⟩ because of opposite signs. Case (iii) leads to a large contribution to the sign problem by |D_i⟩ even though it might be unimportant in the FCI expansion. This indicates that compact wavefunctions due to the orbital representation can on the same note cause a more severe sign problem.

This highlights the general challenge we face, namely, to successfully overcome the sign problem, we may need to abandon the notion of seeking orbital representations that lead to maximally compact wavefunctions.

ΔE^stoq can be modified by changing the basis of the many-particle Hamiltonian. However, finding and applying these many-particle basis rotations that could potentially reduce the gap, in general, are equally hard or even harder than overcoming the sign problem itself.^39,40 Single-particle basis rotations can reduce the sign problem, but this can be used only in very limited cases.⁴¹ 

For the FCIQMC method, however, we find empirically that ΔE^stoq is not the only defining factor on how difficult it is to overcome the sign problem. Since the number of walkers ultimately determines the necessary computational resources for an FCIQMC simulation, the minimum number of walkers N_min to obtain the sign-coherent solution is the important algorithmic quantity. We also define the ratio $s=Nmin/|H|$ as the FCIQMC-related strength of the sign problem. Examples of the relationship between s and N_min with respect to ΔE^stoq are shown in Fig. 4.

We find that s is a monotonic function of ΔE^stoq when moving from half-filling to one hole for a given system size. Increasing the length of the ladder will, however, increase ΔE^stoq but decrease s. On the other hand, when looking purely at N_min, we find that it behaves exponentially with respect to ΔE^stoq for a given filling. The one-hole systems, however, systematically have a higher E^stoq even though a larger but half-filled system might have a larger N_min.

Therefore, we conclude that the exact shape of the wavefunction also influences how computationally expensive it is to overcome the FCIQMC sign problem. Thus, we investigate how guiding wavefunctions that leave ΔE^stoq unchanged influence the strength of the FCIQMC sign problem.

Previously, we have seen that the FCIQMC-related strength of the sign problem s does not only depend on ΔE^stoq but is also system-dependent (see Fig. 3). Hence, instead of modifying ΔE^stoq, we can also make the wavefunction more compact by applying a simple diagonal guiding wavefunction |Ψ_g⟩. Here, we will use a Gutzwiller-like wavefunction ansatz.

It is common practice in various QMC methods to use a guiding wavefunction to improve convergence and reduce stochastic fluctuations.^42–45 We define the Gutzwiller wavefunction as

where |Ψ_HF⟩ is the Hartree–Fock determinant.⁴⁶ |Ψ_HF⟩ spans all determinants in the real-space basis, making the evaluation of ⟨D_i|Ψ_Gutzwiller⟩ require additional computational overhead. To this end, we simplify |Ψ_Gutzwiller⟩ to

where |D_i⟩ are the real-space Slater determinants. The overlaps with another real-space determinant |D_j⟩ are, therefore, simply given by

In FCIQMC, it is applied by scaling the spawning attempts from determinant |D_i⟩ to |D_j⟩ by

In contrast to using full |Ψ_Gutzwiller⟩, these weights are cheap to calculate as the diagonal elements are required anyway for the death/cloning step.

The matrix elements of the matrix that is sampled when weighing the spawns like that are

This is a similarity-transformed version of the original Hamiltonian; thus, the spectrum is conserved. Since the transformation is diagonal and the matrix elements $H̃ij$ are merely a scaled version of H_ji, it is true that

so $Ĥ̃stoq$ is also a similarity-transformed version of $Ĥstoq$⁠. Therefore, the spectrum of $Ĥstoq$ is conserved, too. Thus, ΔE^stoq is left unchanged in the importance-sampled Hamiltonian.

On the other hand, the distribution of the weights in the ground-state wavefunction expansion that is sampled in an FCIQMC run is changed. The ground-state wavefunction of $Ĥ̃$ is given by

where C_i are the real-space wavefunction coefficients of the untransformed wavefunction. For positive values of g, the transformation increases the compactness of the wavefunction, i.e., the ℓ₁ norm of the wavefunction is concentrated on fewer determinants since determinants with larger diagonal elements (i.e., a larger number of double occupancies) are suppressed. The effect on a wavefunction is shown for a small 2 × 3 Hubbard ladder in Fig. 5 where exact diagonalization is trivial.

When running an FCIQMC simulation with a walker number N_w smaller than N_min, the shift S is a biased energy estimator. Since for N_w < N_min the sign problem is unresolved, the sampled FCIQMC wavefunction |Ψ(τ)⟩ is some superposition of the equally valid solutions |Ψ⟩ and −|Ψ⟩. Therefore, the population on the reference determinant |D₀⟩ also does not have a constant sign. Thus, E_proj on average does not converge to the correct ground-state energy and cannot be used as an energy estimator.

The shift S is not only used to control the population but can also be used as an energy estimator. For a large enough population N_w ≥ N_min, the average shift $S̄$ converges to the exact ground state energy E₀. For N_w < N_min, the Hilbert space is not occupied by enough walkers simultaneously such that enough annihilations can take place because opposite sign walkers do not meet in the same iteration on the respective determinant. A sign-incoherent wavefunction is sampled in this case. Because of only partially occurring annihilations, the Hilbert space is overpopulated with respect to the fermionic ground-state wavefunction |Ψ₀⟩ but underpopulated with respect to the stoquastized version $|Ψ0stoq〉$ as some but not all necessary annihilations take place on average. Therefore, the average shift, which is sensitive to the number of annihilated walkers as it is the population control parameter, is in between the two respective ground-state energies in these cases, $E0stoq<S̄<E0$⁠.

Applying a guiding wavefunction such as |Ψ_g⟩ that compactifies the ground-state solution and resembles the true solution |Ψ₀⟩ can greatly improve the convergence of $S̄$ with respect to N_w and, thus, lower N_min. It increases the number of annihilations and in doing so enables the algorithm to resolve the sign problem for a lower number of walkers. Figure 6 shows the convergence curve for a 2 × 8 ladder system at half-filling and with one hole for different g. It also displays the mean number of annihilations for the half-filled system as a function of N_w. Also shown is N_min as determined using the fixed-N₀ scheme with N₀ = 50. As expected, it is situated at a walker number slightly above where convergence is reached.

It is a natural question to ask what the effect of an approximate guiding wavefunction is compared to that using the exact ground-state wavefunction |Ψ₀⟩ if known. Figure 7 shows the dynamics of the walker population and the averaged values for N_min in a fixed-N₀ run for a small 2 × 3 Hubbard system. We compare runs where |Ψ_g⟩ with increasing g and |Ψ₀⟩ are applied as guiding wavefunctions. Increasing g only ever decreases the plateau N_min down to a minimum that is reached when |Ψ₀⟩ is applied. Using |Ψ_g⟩ that is only an approximation to |Ψ₀⟩, however, introduces large fluctuations, which increases the time needed to reduce statistical errors in the energy estimators to the desired level and eventually makes the algorithm unstable, especially when aiming at low N₀ in the fixed-N₀ mode.

We investigate the effect of the Gutzwiller-like guiding wavefunction on the stochastic noise further, this time using the larger 2 × 8 ladder system. As shown in the top panel of Fig. 8, increasing g only ever decreases N_min. Furthermore, we calculate the standard deviation of the shift energy estimator for FCIQMC calculations conducted at their respective N_min for g values ranging from 0.0 (no importance sampling) to 0.2 (see Fig. 8, bottom panel).

For equal CPU time—which corresponds to proportionately more iterations in imaginary time for lower N_w—the lowest standard deviation σ is reached with importance sampling at around g_opt ≈ 0.1. For larger g, σ sharply increases, making exact calculations unfeasible at some point. This effectively puts a boundary on the applicability of |Ψ_g⟩ when deviating from guiding wavefunctions that resemble the exact solution. In some cases, it can be useful, however, to exceed g_opt when one can accept larger statistical deviations in order to further lower N_min to make a calculation affordable.

We can now use the observation that the sign problem in 1D systems is not size-extensive and the possibility to reduce N_min in problems with a size-extensive but weak sign problem, such as the 2D ladder systems (both at half-filling and with one hole). We can use this to calculate the fundamental many-particle gaps in these systems. The fundamental gaps of the Hubbard model are given by

where E(hf) is the energy of the half-filled system and E(1e/h) is the energy with one excess electron or one hole, respectively. Since for the Hubbard model, E(1e) = E(1h) + U, it is only necessary to calculate the half-filled and the one-hole system.

In Fig. 9, we show ΔE^fund for chains up to 102 sites in both open and periodic boundary conditions using the NECI code.⁴⁷ The calculations involved in open boundary conditions are all sign-problem-free. The calculations for the periodic boundary case involve sign problematic calculations; however, these have a non-size-extensive sign problem. Thus, calculations of more than 100 sites are doable. Figure 10 shows the ladder systems with 2 × 6, 2 × 8, 2 × 10, and 2 × 12 sites for both U/t = 8 and 16 calculated with FCIQMC. In the ladder systems, the sign problem is size-extensive. Therefore, the number of walkers to overcome the sign problem scales exponentially with system size. Importance sampling, however, allows for the convergence of the 2 × 12 with an affordable number of walkers. The results and the minimum walker numbers for U/t = 8 are given in Table I. In the most challenging system, the 2 × 12 ladder system with one hole at U/t = 8, with g = 0.25, we were able to reduce N_min from at least several billion walkers—which well exceeds the applicability of our available hardware—to roughly 300 × 10⁶ for which we can calculate the ground-state energy on 32 nodes with 640 cores in 96 h with an error on the order of 10⁻⁴.

In Fig. 11, we compare the FCIQMC result–which is numerically unbiased within stochastic error bars when accounting for the (very small) population control bias—for the 2 × 12 system at U/t = 8 with one hole with results from the density-matrix renormalization group (DMRG).⁴⁸ The DMRG numbers are shown as a function of 1/M, where M is the bond dimension of the matrix product state (MPS). Fiedler-type ordering of the sites is used. The MPS representation is known to be optimal in the 1D case because of the low entanglement of the ground state. In systems that resemble 1D systems, such as the ladder systems, we have been dealing with in this paper, this advantage still carries over. We see good agreement between the FCIQMC result for 3 × 10⁸ walkers and g = 0.25 and DMRG for M ≥ 4000.

In this paper, we have described a certain class of Hubbard systems that show very weak sign problems. Thus, they can be treated in FCIQMC numerically exactly without the initiator approximation. In the special case of 1D systems, the sign problem is non-size-extensive, i.e., systems at various fillings can be treated for system sizes of up to 102 sites when accounting for the population control bias. 2 × ℓ ladder systems show a size-extensive yet very weak sign problem. Using importance sampling with an easy-to-evaluate Gutzwiller-like guiding wavefunction, however, can reduce the required number of walkers such that systems well beyond the scope of exact diagonalization can be treated numerically exactly within stochastic error bars.

In an upcoming publication, we will discuss controlled approximations in systems beyond what is possible in a numerically exact manner. Importance-sampled FCIQMC will be used as a basis to treat very-weak sign problem subspaces of larger systems. Furthermore, using importance sampling with more sophisticated guiding wavefunctions, if evaluated efficiently, could potentially reduce the minimum number of walkers further, expanding the scope of treatable systems.

The authors would like to thank Michael Willatt for fruitful discussions. The authors acknowledge the Max Planck Society for support.

The authors have no conflicts to disclose.

Niklas Liebermann: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); writing – original draft (lead); writing – review & editing (lead). Khaldoon Ghanem: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Software (supporting). Ali Alavi: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (supporting); writing – original draft (supporting); writing – review & editing (supporting).

All FCIQMC calculations were conducted using NECI (https://github.com/ghb24/NECI_STABLE). DMRG calculations were conducted using BLOCK (https://github.com/sanshar/Block). The data that support the findings of this study are available from the corresponding author upon reasonable request.

The 1D Hubbard model is sign-problem-free in the following cases:

1. for open boundary conditions and

2. for periodic boundary conditions for odd n_↑ and n_↓.

This can be easily seen in the usual fermionic second quantization formalism. To this end and without loss of generality,⁴⁹ we will write a Slater determinant |D_i⟩ with the fermionic creation operators $ĉp,σ†$ at site p acting onto the vacuum state |⟩ in the convention of spin-first ordering as

with a positive sign. n_↑ (n_↓) are the number of up-spin (down-spin) electrons. α_p and β_q number the spatial sites in ascending order. The Hubbard Hamiltonian only allows for hoppings of the same spins to nearest-neighboring sites, so without loss of generality, an excitation of an ↑-electron from site j to a site j ± 1 is given by

(where j ± 1 is not occupied with an ↑-electron already). Since both $ĉj±1†$ and $ĉj,↑$ were commuted with the same number of $ĉ↑†$ operators, trivially, there is no sign change. For an excitation of a ↓-electron, the same applies. This proves that all 1D Hubbard systems with open boundary conditions are sign-problem-free. With an excitation that exploits periodic boundary conditions, however, we get

$ĉℓ,↑$ had to be commuted with n_↑ − 1 creation operators, whereas $ĉ1,↑†$ did not have to be commuted at all. Thus, there is no sign change only for odd n_↑. For a ↓-electron, we get

In addition, here, n_↓ has to be odd for no sign change to occur.

Without loss of generality, let us look at two determinants |D_i⟩ and |D_j⟩ in a system with an even number of ↑-electrons and an odd number of ↓-electrons. Between |D_i⟩ and |D_j⟩, a single ↑-electron is moved from site i to site j,

The two main contributing pathways consist of one moving ↑-electrons not exploiting periodic boundary conditions and one moving them exploiting the boundary conditions. In the FCIQMC language, according to Eq. (3), a walker sitting on |D_i⟩ is creating

walkers on |D_j⟩, where n_docc is the number of doubly occupied sites shared by |D_i⟩ and |D_j⟩, |j − i| is the number of sites between site i and site j, and $n↓i,j$ is the number of ↓-electrons between i and j. In contrast, exploiting periodic boundary conditions leads to

On the one hand, the contribution of these pathways to the fermionic solution is given by P_non-periodic + P_periodic. On the other hand, the contribution to the stoquastic solution is $Pnon-periodic+Pperiodic$⁠. Therefore, we can conclude that the stoquastized gap due to these contributions is largest when P_non-periodic = −P_periodic because of maximal cancellation. One deviates from the equality if (i) n_sites is large, (ii) $n↓−n↓i,j$ is large, and (iii) U is large, implying a reduction of the sign-problematic contribution.