Sign-free stochastic mean-field approach to strongly correlated phases of ultracold fermions

We propose a new projector quantum Monte-Carlo method to investigate the ground state of ultracold fermionic atoms modeled by a lattice Hamiltonian with on-site interaction. The many-body state is reconstructed from Slater determinants that randomly evolve in imaginary-time according to a stochastic mean-field motion. The dynamics prohibits the crossing of the exact nodal surface and no sign problem occurs in the Monte-Carlo estimate of observables. The method is applied to calculate ground-state energies and correlation functions of the repulsive two-dimensional Hubbard model. Numerical results for the unitary Fermi gas validate simulations with nodal constraints.

2. The sign-free stochastic mean-field scheme. For a system of fermions interacting through a binary potential, we first introduce a set of hermitian one-body operators s Â ( )  + are the Fermi creation and annihilation operator in a single-particle mode i . Let us now consider two N-particle Slater determinants ψ , ϕ of orbitals satisfying the biorthogonality relations Appendix that the Slater determinant ϕ transforms according to: A would represent physically one-particle-one-hole and two-particle-two-holes excitations. After an infinitesimal time step dτ , during which ϕ evolves to , the last term of the expansion (2) causes departures of the propagated state from a single Slater determinant. However, introducing random fields according to the Hubbard-Stratonovich transformation can linearize the dynamics [19][20][21]. The detailed derivation is reported in Appendix. Finally, the exact imaginary-time evolution is recovered through the stochastic average of Brownian trajectories in the subspace of Slater determinants states ϕ : is the average over a random functional; the evolution of S and the orbitals n ϕ is governed by the following equations in the Itô sense: The s dW are infinitesimal increments of independent Wiener processes: τ . Its overlap with the exact many-body ground-state g Ψ can be obtained from the representation (3) in the limit of large τ , provided that the trial state ψ is not orthogonal to g Ψ : where g E is the ground-state energy and where the phase of ψ can be chosen to have and no walker can cross the exact nodal surface if S is real. In contrast, the standard auxiliary-field projector quantum Monte-Carlo method [22], as well as our previous stochastic mean-field scheme [16,17], generally leads to Slater determinants ϕ whose overlap ϕ ψ exhibits a varying sign. As a consequence, a walker can reach the nodal surface and it then generates stochastic paths that do not contribute to the ground-state: as long as Such trajectories only increase the statistical error and are responsible of the sign-problem. Thus, one is forced to perform the constrained-path approximation [23] where walkers are eliminated "by hand" when their overlap with a ground-state ansatz wave-function becomes negative.
Numerically, the stochastic differential equations (5) are solved, in the Stratonovich form, by an embedded Runge-Kutta (5,4) algorithm with adaptive stepsize control. We take ( ) ψ ϕ = 0 as initial condition and the spreading of the weights ( ) S exp is avoided through standard population control techniques. In practice, we use the stochastic reconfiguration method [24] that deals with a fixed-number of walkers ϕ among which some are killed and others are cloned according to their relative weight in the population. Observables are estimated from the representation (3) of the manybody state. For instance, the ground-state energy E g is be obtained according to where the amplitudes are now all real positive. Moreover, when the Slater determinant ansatz ψ and the ground-state share a common symmetry, the stochastic paths are automatically projected onto this symmetry sector in the estimate (7). Otherwise, the sampling can be improved by projection techniques [25,26]. These conclusions also hold true for any observable commuting with the Hamiltonian. In other cases, one obtains an approximate ground-state expectation value, known as the mixed estimate [27]. It can be corrected by the following extrapolated estimate that is one order better in the difference g Ψ − ψ [27]: Note that these observable estimates are biased in standard projector quantum Monte-Carlo methods by the nodal constraints introduced to circumvent sign problems. The stochastic Hartree-Fock approach (3)(4)(5) removes this systematic error. The exact expectation value would require obtaining the ground-state density matrix by a double propagation in imaginary-time: For such calculations, our scheme can be extended by expanding the density-matrix in terms of dyadics that are formed by biorthogonal Slater determinants, both undergoing a Brownian motion similar to Eq. (5). We emphasize that the method then becomes the equivalent, at fixed particle-number, of the recent Gaussian Monte-Carlo technique [25,28], which is more suited for thermodynamical studies in the grand-canonical ensemble.

Application to ultracold atomic
Here + σ . In other cases, one experiences severe sign problems that practically prohibit studying large lattices, strongly interacting systems or open shells configurations [22]. In contrast, our new stochastic Hartree-Fock scheme (3-5) does not manifest explicit sign problems regardless of the lattice topology, band filling and sign of the interaction. Indeed, a quadratic form (1) can be recovered from the Hamiltonian (9) by using the local density or magnetization depending on the sign of the interaction parameter U :   Fig. 1 for two-dimensional lattices with a hole doping 125 . 0 = δ from half-filling. This density generates the most important sign problem in auxiliary-field approaches [22]. The antiferromagnetic mean-field state is very efficient and the exact ground-state energy of the 4 4 × lattice at t U 4 = is recovered to within less than 0.2% if the sampling is improved by projecting onto the spin-singlet sector. For the larger on-site repulsion t U 8 = , the stochastic Hartree-Fock estimate of the ground-state energy is The simulations on the 8 8 × lattice, presented on the right panel of Fig. 1, are in agreement to less than 0.5% with the constrained-path approach [23]. This is fully consistent with the error usually observed in constrained-path calculations of the ground-state energy on small clusters [23]. But, we emphasize that such small discrepancies can originate from the numerical error in the integration of the stochastic differential equations (5) and from the unavoidable bias that is introduced by the population control algorithm. In the half-filling limit, our calculations confirm the emergence of an antiferromagnetic phase [31], as shown in Fig. 2    We finally address the unitary Fermi gas limit. In this ideal regime of strong interaction via a two-body potential of zero range and infinite scattering length, fermions are among the most intriguing physical systems since they are believed to exhibit universal many-body states. For instance, at zero temperature, the energy must be a universal fraction η of the Fermi energy that is the only relevant energy scale in the system. From experiments with trapped atomic gases [5][6][7], the measured values for this ratio η vary from 10 Here periodic boundary conditions are assumed in each direction; are the matrix elements of the single-particle kinetic energy operator in the representation position; M is the atomic mass, l denotes the grid spacing and K 44275 .
is a numerical constant. In the unitary limit, a goes to infinity but the coupling constant on the lattice remains finite and negative, so that the gas clearly experiences attraction. Our sign-free simulation method with the model Hamiltonian (12), transformed as in Eq. (10), has been checked from the known solutions of the two-and three-body problem in an harmonic trap at the unitarity point [35,36]: in all cases, the discrepancy does not exceed one percent. For different systems, up to 42 = N atoms on a 8 8 8 × × lattice, we plot in Fig. 3     We acknowledge fruitful discussions with Y. Castin, R. Frésard and F. Gulminelli.

Appendix. Derivation of the sign-free stochastic Hartree-Fock equations.
Consider, at a given imaginary-time (omitted for simplicity), two N -particle Slater determinants ψ , ϕ with biorthogonal orbitals:.
where Ω is a non-singular matrix and k J a k k d d × bi-diagonal matrix: With the help of (A. 7), only values of the integer l ranging from 1 to N give a non-zero contribution, and for N k ≤ the single-particle wave-function k k ϕ ω = belongs to the Fermi sea ϕ , so that ϕ δ ϕ ω ω kl l k a a = +~. One is finally left with: with the following variation of the Hartree-Fock orbitals: Here, we have used the idempotency of R that implies ( ) 0 . Note that he evolution process (A.17) preserves the biorthogonality constraints (A.1). Therefore, the propagation scheme (A.16) can be iterated for an arbitrary number of imaginary-time steps τ d and one is naturally led to the stochastic mean-field approach (3)(4)(5). The results have been averaged over several hundred of runs of 100 trajectories. When not shown, statistical error bars are smaller than the symbol size. The black line gives on the left the exact groundstate energy [32] and on the right the constrained-path Monte-Carlo result [23].  is the ratio between the mean-energy of at "time" τ and the ground-state energy of the non-interacting gas on the lattice. The trial state is a spin-singlet Slater determinant for the free gas. We show the average result over many runs of 100 paths. We use a 6 6 6 × × lattice, except for N ≥ 40 where the calculations were performed on a 8 8 8 × × grid. The imaginary-time τ is expressed in units of 2 2 h ml , where l is the lattice spacing.