Few-body physics on a space-time lattice in the worldline approach

We formulate the physics of two species of non-relativistic hard-core bosons with attractive or repulsive delta function interactions on a space-time lattice in the worldline approach. We show that worm algorithms can efficiently sample the worldline configurations in any fixed particle-number sector if the chemical potential is tuned carefully. Since fermions can be treated as hard-core bosons up to a permutation sign, we apply this approach to study non-relativistic fermions. The fermion permutation sign is an observable in this approach and can be used to extract energies in each particle-number sector. In one dimension, non-relativistic fermions can only permute across boundaries, and so our approach does not suffer from sign problems in many cases, unlike the auxiliary field method. Using our approach, we discover limitations of the recently proposed complex Langevin calculations in one spatial dimension for some parameter regimes. In higher dimensions, our method suffers from the usual fermion sign problem. Here we provide evidence that it may be possible to alleviate this problem for few-body physics


I. INTRODUCTION
Computing the properties of quantum systems containing fermions remains challenging especially when perturbative techniques begin to fail. Even in the case of few-body physics, where each particle is described by a large dimensional vector space, the free and the interacting parts of the Hamiltonian may be diagonalized by two very different basis vectors and the ground state in a given particle-number sector may be severely entangled in both these bases with no apparent small parameters. This problem is even more severe in quantum field theories, where these particles are created out of a vacuum that can itself be non-trivial, like in Quantum Chromodynamics (QCD). For this reason, studying the properties of low-mass hadrons remains a daunting challenge in lattice QCD [1].
In the context of QCD, an alternative approach has become exciting in recent years and is based on ideas of a low-energy effective field theory formulated using the symmetries of QCD and the fact that the vacuum breaks the chiral symmetry spontaneously [2,3]. This chiral effective field theory (χ-EFT) is constructed using nucleons and pions as the low-energy degrees of freedom. The interactions are described by local operators constructed from the nucleon and pion fields, organized as a power series in the ratio of relevant energy scales, which acts as the small parameter. But even at leading orders in the effective field theory, computing the properties of low lying hadrons and nuclei within χ-EFT can require non-perturbative calculations. While analytic techniques based on resummations of Feynman diagrams are useful for up to three particles [4], numerical approaches especially based on Monte Carlo methods become necessary for higher number of particles [5,6].
Non-perturbative state of the art Monte Carlo methods for few-body problems are based on auxiliary field techniques and fixed node approximations. Many variants of the algorithms have been developed over the years and for further details we refer the reader to recent reviews of the subject [7][8][9]. While these auxiliary field quantum Monte Carlo (AFQMC) methods have a clear advantage in certain parameter regimes, they also exhibit several limitations in other regimes [10]. Difficulties of the method become apparent in one spatial dimension, which has become experimentally interesting in recent years, thanks to our ability to design and control ultracold quantum gases confined to optical traps [11,12]. Also, many interesting quantum phenomena in higher dimensions have analogues in one spatial dimension [13][14][15][16][17]. Motivated by this, the AFQMC method was recently used to study two species of fermions in one spatial dimension interacting through a delta function interaction [18][19][20].
One major limitation of the auxiliary field approach is that it suffers from sign problems in the presence of repulsive interactions or mass-and spin-imbalanced systems away from half filling. This is true even in one spatial dimension. In order to explore a solution to this problem, the auxiliary field technique was combined with the complex Langevin (CL) method to overcome the sign problem in Ref. [21]. Recently, this approach was also extended to higher dimensions [22]. Unfortunately, it is well known that CL methods may have uncontrolled systematic errors and can converge to wrong results [23]. In fact, the authors of Ref. [21] suggested caution for their results, especially in the repulsive case where the CL approach showed fat-tailed distributions. However, the authors wondered if the flattening of the ground-state energy as a function of the strength of repulsion, observed using the CL approach, was a sign of some interesting non-perturbative physics. The only way to be sure is to compare the results with other reliable methods.
One spatial dimension is an excellent place to test methods like CL, since entanglement is greatly reduced and a variety of other methods that do not suffer from sign problems are usually available. For example, exact analytic calculations based on the Bethe ansatz are possible in some special cases [24][25][26][27]. With open boundary conditions or odd number of fermions with periodic boundary conditions, sign problems are absent in the worldline formulation [28,29]. Recently the two-dimensional lat-tice Thirring model with open boundary conditions was formulated using the worldline method [30] to provide benchmark results to test the Lefschetz thimble approach [31]. Thus, we should also be able to test the recent CL results of Ref. [21] using a similar worldline approach.
Motivated by this, we construct a general worldline approach to study quantum mechanics of hard-core bosons in any dimension. We show that we can use worm algorithms to update the worldline configurations efficiently in any particle-number sector by tuning the chemical potential carefully [32]. While worldline methods to study bosonic quantum field theories are well known by now [33][34][35], and have been used in several studies so far [36][37][38][39][40][41], their applicability to bosonic quantum many-body physics has remained relatively unexplored [42,43]. Thus, our work can be viewed as an attempt to fill this gap. Our method can easily be extended to fermions if we can compute the fermion permutation sign accurately. In one spatial dimension, fermions are identical to hard-core bosons up to boundary effects. Hence, our method can be used to check the recent results of Ref. [21]. We find that the CL method yields incorrect results as the repulsive coupling strength grows, implying that the observed flattening is unphysical and an artifact of the method.
We can easily adapt our approach to study fermionic particles even in higher dimensions by treating the fermionic permutation sign as an observable, but it becomes difficult to compute it accurately. As expected, this observable suffers from a severe signal-to-noise ratio problem, especially at low temperatures and when the number of particles becomes large. However, we provide evidence that at intermediate temperatures and with a small number of particles (N ∼ 10) we may be able to beat the signal-to-noise ratio. This may allow us to explore new and interesting questions in few-body physics that are difficult to answer with the AFQMC method.
Our paper is organized as follows. In Section II we discuss the lattice Hamiltonian formulation of two species of hard-core bosons with a contact interaction. In Section III, we construct the worldline formulation of the problem on a space-time lattice and provide details of our worm algorithm in Section IV. In Section V, we discuss how we can study fermions by measuring the fermionic permutation sign. We also provide evidence that the fermion sign problem is mild in three dimensions for up to ten particles at an intermediate temperature. In Section VI, we present our results for fermions in one spatial dimension and consider two cases: the mass-balanced case and the mass-imbalanced case. In the mass-balanced case, we show that our results are in agreement with the exact results obtained using the Bethe ansatz. We also show that in both cases the flattening of the ground state-energy observed in the CL calculations is absent in our approach. In Section VII, we discuss a limitation of the traditional approach used to extract the ground state energy and suggest a complementary method that is more reliable. In Section VIII, we provide evidence that we can also compute fermionic ground state energies in 3 + 1 dimensions by reproducing a benchmark calculation performed several years ago in Ref. [44] and present our conclusions in Section IX.

II. LATTICE MODEL
The physics of two species of non-relativistic hard-core bosons that we consider in this work can be represented through the lattice Hamiltonian where c † i,σ , c i,σ create and annihilate bosons of species, say, spin-up (σ = ↑) or spin-down (σ = ↓), at the site i on a d-dimensional spatial lattice andα represents unit vectors in the d spatial directions. The site occupationnumber operators are defined as N i,σ = c † i,σ c i,σ . The parameters t σ = 1/(2m σ a 2 ) give the hopping strength of the particles in terms of the masses of the particles m σ and the lattice spacing a. The parameter U is the bare interaction strength. In addition, we impose the hard-core boson constraint on the states of the Hilbert space, which means each site can either be empty or contain a single boson of a particular species. In this work, we will study finite spatial boxes of size L X with periodic boundary conditions, so that the physical box size is given by L = L X a. If the particles are taken to be fermions instead of hard-core bosons, the lattice model (1) is known as the Hubbard model.
The lattice model (1) has a naive continuum limit as a → 0. In d = 1, this limiting procedure leads to the continuum theory with a delta function interaction where g = U . However, in higher dimensions (d > 1) the problem of the continuum limit is more subtle and the framework of effective field theory becomes necessary to implement it. For example, with equal-mass fermions in three dimensions, the coupling U (a) can be tuned to the unitary fixed point to get an interacting field theory in the continuum limit. With bosons or additional quantum numbers, we can get Efimov physics which necessitates additional counterterms to renormalize the theory [4]. For this reason, we confine the discussion of the continuum limit to one dimension. In higher dimensions we will simply view (1) as a lattice model and set a = 1.
In this work, we will compute the ground-state energy E 0 N ↑ ,N ↓ , which is the lowest eigenvalue of the Hamiltonian in Eq. (1) in the subspace with particle numbers N ↑ and N ↓ . One way to accomplish this is to simply compute the average energy at very low temperature (β → ∞), where we use the modified lattice Hamiltonian, with chemical potentials µ σ for the two particle species, and the partition function to define the expectation value. If we compute the trace in a fixed particle-number subspace, the chemical potential terms drop out and we indeed get E = E 0 N ↑ ,N ↓ in the limit β → ∞. Thus, we need a method to compute E reliably for large values of β.
In the next two sections, we will find an expression for E in the worldline formulation and construct a worm algorithm to compute it efficiently. Worm algorithms work by adding and removing particles which then updates the worldline configuration. If energetics do not favor this, algorithms can develop exponentially long autocorrelation times. Hence, for efficient sampling in a fixed particle-number sector it is important to tune the chemical potentials carefully. To understand why this is the case, let us consider the energy of the ground state containing N ↑ and N ↓ particles in the presence of a chemical potential, which is given by If we can tune µ σ to critical values µ c σ such that then worm algorithms can sample worldline configurations in the particle-number sectors (N ↑ , N ↓ ) and (N ↑ + 1, N ↓ + 1) very efficiently even at very large values of β. We can monitor this by computing the average particle number, and making sure that it fluctuates between the sectors (N ↑ , N ↓ ) and (N ↑ + 1, N ↓ + 1) even when β is large [36]. Such fluctuations are crucial to the efficiency of our algorithm.

III. WORLDLINE FORMULATION
Let us now construct the worldline formulation of the problem. We first write the Hamiltonian as H µ = H d −H h , a sum of a diagonal term and a hopping term, where We then expand the partition function as which can be viewed as a hopping parameter expansion in continuous time. Since we do not truncate the sum over k in Eq. (11), there is no approximation involved. Such an approach to write partition functions in continuous time is well known for hard-core bosons [45] and fermions [46]. However, to develop the worm algorithm, it is convenient to discretize the time integrals by dividing β into L T imaginary time steps of width ε, such that β = εL T . If we then compute the trace in the occupation number basis, we can approximate the partition function as a sum over worldline configurations C of both species of particles on a space-time lattice. We write this as where Ω(C) is the Boltzmann weight of each worldline configuration. Figure 1 gives an illustration of C on a 1 + 1 dimensional space-time lattice. In the next section, 1. An illustration of a worldline configuration C with N ↑ = 1 and N ↓ = 1. The dots represent space-time lattice sites and the bold solid lines show the worldlines of the two particles. The interaction between the particles is shown as a wiggly temporal bond. The Boltzmann weights associated with lattice sites containing particles are also shown.
we construct a worm algorithm to update such worldline configurations. We usually perform calculations at several values of ε and then extrapolate to the continuous time limit. We always find that the time-discretization errors are linear in ε at leading order (see Figs. 13 and 14). In principle, this extra work can be avoided by directly taking the continuous time limit of the worm algorithm itself [45]. The Boltzmann weights Ω(C) can be obtained as a product of local weights associated to each space-time lattice site, where W i,t (C) is obtained from the local worldline configuration C at the space-time site (i, t) and is a product of either W σ e , W σ h or W σ m for each species and W I that takes into account particle interactions. The allowed configurations at a site for each particle species are shown in Fig. 2. Either there is no particle or one particle that comes into the site and leaves it. If there is no particle, the site weight for that species is W σ e = 1. On the other hand, when there is a particle, we choose the weight to depend only on its outgoing direction. If the particle hops to the neighboring spatial site, the weight is If, instead, the particle moves forward in time, the weight is The interaction among the particle species is taken into account through the weight if both particles hop forward in time together. Otherwise we set W I = 1. With these definitions, we can express the weight of a configuration as where N σ is the number of σ-particles in the configuration C, n σ h is the number of hops and n I is the number of interacting temporal bonds.
It is possible to compute a variety of observables easily in the worldline formulation. For example, the average energy defined in Eq. (3) can be expressed in the worldline formulation, up to O(ε) errors, as where E(C) is the energy of a worldline configuration C, defined as The allowed worldline configurations C on a 1+1 dimensional space-time lattice site for one particle species.
Using the weights defined in Eqs. (14) to (16), the weights of the configurations from left to right are (top row) 1, Wm, Wm, In addition to these weights if both layers have the weight Wm we multiply the product with WI .
This expression for E(C) can be derived from Eq. (17) by taking the appropriate derivative with respect to β.
The average particle number for each species defined in Eq. (8) is straightforward since worldline configurations have a well-defined particle number N σ (C). Thus, we get We can also measure ratios of partition functions, like the one we define later in Eq. (35), by designing an appropriate sampling method during the Monte Carlo update.

IV. THE WORM ALGORITHM
It is possible to develop worm algorithms to update the worldline configurations C of the type shown in Fig. 1 [47,48]. During the worm update, we pick each particle species and update its worldline configuration while keeping the worldline of the other species fixed. To perform the update, we pick a space-time site at random and create a defect in the worldline configuration in the form of either a particle creation or annihilation event. We choose the worm head as the position of the annihilation operator and the tail as the position of the creation operator. We then move the head on the space-time lattice, while keeping the location of the tail fixed, using local moves that obey detailed balance. The worldline configuration gets updated on each site the worm visits. The update ends when the worm head meets the tail and removes the defect. The particle number can be monitored during the worm update and local rules can be chosen so as to sample particle numbers within a fixed range.
For a given worldline configuration, we define the local configuration at a space-time site as the incoming and the outgoing directions of the particle. In addition we also include the information about how the worm enters or leaves that site. Each box in Figs. 3 and 4 represents such a local configuration with the worm entering the site. Local configurations with the worm leaving the site can be constructed from these by reversing the direction of the worm arrow. We then group configurations that can transform into each other under local rules which satisfy detailed balance.
To understand this procedure better, let us consider local updates that begin or end a worm update, shown in Fig. 3 for a 1 + 1 dimensional lattice. When the worm update begins, the incoming worm direction is shown as a diagonal arrow entering the site. The outgoing worm direction could be along any of the neighboring space-time lattice sites as long as that move is allowed. For detailed balance to work, the worm update should also be allowed to end through the reverse process. There are two classes of such begin-end updates based on whether the first site contains a particle or not. If the site contains a particle (Fig. 3a), then a creation operator is introduced on the site (i.e., the site becomes the worm tail) and the worm head that contains the annihilation operator moves to the neighboring site through which the particle came to the first sit (i.e., with probability one the outgoing direction is chosen to be the backward direction the worldline). For detailed balance to be satisfied the reverse process is also chosen with probability one (i.e., when the worm head enters the site that contains the tail, the worm update ends). Thus, these two local pair of configurations are grouped together. There are seven possible such pairs of begin-end updates in 1 + 1 dimensions, as shown in Fig. 3a. In higher dimensions there would be more.
It is also possible that the worm update begins at an empty site (Fig. 3b). In such a case, a new particle is proposed to be created on that site (i.e, a creation operator is placed on the site which becomes the worm tail) and the worm head moves either to one of the neighboring spatial sites or upwards to the neighboring temporal site. In 1 + 1 dimensions, the empty site can thus transform into three possible local configurations. These four configurations are grouped together and shown in Fig. 3b. We assign probabilities for moves within this set of four local configurations such that detailed balance is satisfied. Details on how these probabilities are chosen are discussed below.
Between the begin and end updates, the worm head moves through space-time lattice sites. Such moves can be classified into three classes, as shown in Fig. 4. The simplest class involves a move where the worm enters the site that already has a particle on it. In this case, the entering direction of the worm cannot be the same as the incoming particle direction. Thus, with probability one, the outgoing worm direction can be exchanged with the incoming particle direction. Eight such pairs of configurations exist in 1 + 1 dimensions and are shown in Fig. 4a. The next class involves a worm moving into an empty site from a spatial neighbor. Then the worm can exit through the forward time direction or through a different spatial direction. This groups 2d + 1 local configurations together in d spatial dimensions. The two possible groups of three configurations in 1 + 1 dimensions are shown in Fig. 4b. The third class involves a worm entering into an empty site along the forward time direction. Then the outward worm direction could be along any of the spatial or forward time directions. In d dimensions, there are 2d + 2 such local configurations that can transform into each other. The only possible such group of four configurations in 1 + 1 dimensions is shown in Fig. 4c.
For the efficiency of the worm method, probabilities for moving the worm head must be chosen so as to avoid a "bounce" as much as possible [49]. A bounce occurs when the local configuration does not change (i.e, the Case I: W2 > W1 + W3 probability to simply reverse the incoming worm direction wins). Note that bounces can be completely eliminated among the pairs of configurations in Figs. 3 and 4, since they have the same weight. In these cases, the local worm update simply toggles the two. However, in the case of groups that contain more than two local configurations, like the four-configuration group in Fig. 3b or the threeand four-configuration groups in Figs. 4b and 4c, we need to make sure the worm bounces are minimized as much as possible. For small values of ε, since two of the weights are very close to one, and spatial hops are rare this is almost always possible.
Let P ab be the transition probability to go from local configurations a to b. To illustrate how we choose P ab to satisfy detailed balance, we consider the example of a three-configuration group in Fig. 4b. The bounce probability is given by The weights of the three configurations are W 1 = 1, W 2 = W m (or W 2 = W m W I ) and W 3 = W h . In this case, W 1 and W 2 are close to one, while W 3 is of the order of ε.
There are three possible orderings of the weights: Case I: Table I, we give the transition probabilities P ab for all these cases.
The above method of constructing transition probabilities is easily extended to other groups of local configurations in our work. For example the group of four configurations in Figs. 3b and 4c have weights W 1 = 1, In these cases we can imagine the configurations with weights W 3 and W 4 as a single configuration with weight W 3 +W 4 = 2W h and again use the probabilities of Table I. The complete set of transition probabilities for this case is given in Table II. We test our algorithm for small lattices in d = 1, 2, 3 dimensions and the results of these tests are given in Appendix A.
where S(C) is the fermion permutation sign of the worldline configuration C. If we define the average fermion permutation sign S as the ratio of the two partition functions in a fixed particle number sector as then the energy observable in the fermionic theory can be defined using the relation where E is the average energy for hard-core bosons computed through Eq. (3). The ground-state energy of the theory with fermions is then obtained in the low temperature limit The situation simplifies in one spatial dimension. It can be shown that with periodic boundary conditions and an odd number of particles of each species, S = 1. In the case of the open boundary conditions or in a trapping potential, the fermions cannot wind around the box and S = 1 again. This implies that the fermionic and bosonic energies are the same, that is, E 0f N ↑ ,N ↓ = E 0 N ↑ ,N ↓ . We use this result in the next two sections when we study one-dimensional problems. In higher dimensions, we expect S < 1. In particular, the average permutation sign will suffer from a severe signal-to-noise ratio problem when β becomes large and in the presence of large number of particles. This is essentially the sign problem coming back to haunt the Monte Carlo method. However, it is interesting to note that our algorithm does does not distinguish the free problem (U = 0) from the interacting cases of attraction (U < 0) or repulsion (U > 0). Our approach also does not differentiate between whether the two species have similar or different masses. The sign problem is severe in all cases, although we find it somewhat milder with attractive as compared to repulsive interactions. In Fig. 5, we plot the average sign S as a function of the particle number N = N ↑ + N ↓ in the mass-balanced repulsive model with m σ = 1 and U = 4 at temperatures such that β(1 − cos(2π/L)) ∼ 1.5. For even N we choose N ↑ = N ↓ and for odd N we use N ↑ = N ↓ + 1. We see that up to N = 10 the sign problem is not very severe at these intermediate temperatures. This gives us confidence that, in combination with ideas of fermion bags [50] and exponential error-reduction techniques [51], we may be able to use this approach to study few-body physics even in higher dimensions. We show some preliminary evidence for this in Section VIII.  [21]. Note that they begin to disagree with exact results around γ 2.5. 2 in the right figure shows that the worm algorithm results are sensitive to the small difference between the lattice and continuum. In the β F → ∞ limit, we recover the exact lattice result as expected.

VI. ONE-DIMENSIONAL SYSTEMS
In this section, we report our results for a wide range of coupling strengths and mass-imbalances in one spatial dimension for spin-balanced systems (N ↑ = N ↓ ). We fix the box size to be L X = 40 and impose periodic boundary conditions. As mentioned in the previous section, an odd number of fermions of each species is equivalent to hard-core bosons and our so method is directly applicable to fermions. We focus on computing the average energy  [21] and begin to disagree with the worm algorithm around γ 2, even in the regime where second-order perturbation theory works well.
E (see Eq. (3)) at a fixed β = 100. In the next section, we discuss how our results change with β for a few cases in more detail and how understanding this dependence is important to extract the ground-state energy accurately. We do perform an extrapolation to zero temporal lattice spacing ε → 0, since these errors can be large. Details of this extrapolation are discussed further in Section VIII and Appendix A.
We denote the total number of particles as N = N ↑ +N ↓ and the particle-number density by n = N/L X . We define the average mass and the mass-imbalance parameter as respectively. We work in units with m σ = 1 and a = 1 and parametrize the interaction strength by γ, which is related to the bare coupling U in the lattice model (1) by To facilitate comparison with literature, we report all energies in units of the corresponding one-dimensional ideal spin-and mass-balanced Fermi-gas ground state energy E FG in the continuum defined as where F = π 2 n 2 /8m is the Fermi energy.

A. Mass-Balanced Systems
The continuum model (2) for fermions with m = 0 is known as the Gaudin-Yang model and can be solved exactly when g > 0 using nested Bethe ansatz [24,25]. This solution can also be extended to the lattice model (1) in one spatial dimension (d = 1) for an arbitrary lattice size L X , when coupling U ≥ 0, and number of particles N = N ↑ + N ↓ as long as m ↑ = m ↓ [26]. Since attractive models on a finite lattice can be related to repulsive models by a particle-hole transformation on one of the spins, we can compute the exact ground-state energies of our lattice model in one spatial dimension when m = 0 for a wide range of the coupling strengths, both attractive and repulsive. This allows us to test our Monte Carlo method for the mass-balanced case against the exact solution. As we show below, we find excellent agreement between the two.
For completeness, let us quickly review the main steps that go into the exact computation of the ground-state energy. We discuss the solution for the lattice model only, since the solution in the continuum can easily be obtained from it. Let us assume we have N ↓ = M spin-down fermions, and N ↑ = N − M spin-up fermions. Let us normalize our energies with t ↑ = t ↓ = t. Let p i = (p 1 , . . . , p N ) and λ α = (λ 1 , . . . , λ M ) be two sets of ascending real numbers. Then the following N + M coupled non-linear equations in the N + M variables {p i , λ α } determine the complete spectrum of this system: (30) where i = 1, . . . , N , α = 1, . . . , M , and I i , J α are specific integers (or half-odd integers) for N, M even (or odd) that uniquely label the energy eigenstate. The energy eigenvalue is then given by E = 2t i (1 − cos(p i )). For the ground state, we must choose I i = − 1 2 (N + 1) + i and J α = − 1 2 (M + 1) + α. The ground state for the attractive Hubbard model can be obtained from the repulsive case using the relation The exact ground-state energy obtained by this procedure for N ↑ = N ↓ = 5 is shown in the left plot of Fig. 6 as a dashed line labeled "Bethe ansatz." We also show results obtained from perturbation theory up to second order. We note that our method, labeled as "Worm algorithm," is able to precisely reproduce the exact results from the Bethe ansatz once we extrapolate to the continuous time limit ε → 0 even at β = 100. For comparison, we also plot the CL results from Ref. [21] on the same plot. We notice that CL reproduces the results quite well for small values of the couplings, but starts to deviate when γ 2.5. To confirm that the worm algorithm reproduces the exact Bethe ansatz calculations in the strong coupling regime, where perturbation theory clearly breaks down, we extend these results to higher values of γ in the right plot of Fig. 6.
There is a small difference between the exact results from Bethe ansatz in the continuum as compared to the lattice, as shown in the left plot of Fig. 7. In the right plot, we show that the worm algorithm is able to resolve this difference. the exact continuum and lattice Bethe ansatz results (at β F = ∞). We also show the worm algorithm results at different values of β F , which indeed converge to the correct lattice answer in the limit β F → ∞.

B. Mass-Imbalanced Systems
In contrast to the mass-balanced case, no exact solution is known for the general mass-imbalanced case (m = 0). However, our algorithm is again very efficient even in such cases. In this section, we compute some spinbalanced results and check our results against secondorder perturbation theory and find good agreement for small values of γ. We show our data using the worm algorithm along with the iHMC results of Ref. [21], unpublished CL data given to us by the authors of Ref. [21], and second order perturbation theory.
The apparent discrepancy at m = 0 between our results and the results from the Bethe ansatz is addressed in Fig. 10.
ground state energy of 5 + 5 particles as a function of the coupling for a high mass-imbalance of m = 0.8. We perform computations at β = 100, which we find to be sufficient to make our point. As can be seen from the left plot of Fig. 8, our results agree very well with second-order perturbation theory for up to γ ∼ 3.0, while the CL results disagree with both perturbation theory and our method when γ 2. In the right plot of Fig. 8, we extend this to the regime of very strong repulsion. It was suggested in Ref. [21] that the flattening of the ground-state energy as a function of γ for strong repulsive couplings, obtained using the CL method, could be a physical effect. However, the disagreement with exact Bethe ansatz calculations at m = 0 (Fig. 6) and with the worm algorithm at a high mass-imbalance of m = 0.8 (Fig. 8) shows that the observed flattening is an artifact of the CL method.
In Fig. 9, we perform a comparison with the results of Ref. [21] for N ↑ = N ↓ = 3 fermions with a fixed attractive coupling γ = −3.0 across a range of mass-imbalances. We observe that imaginary-mass Hybrid Monte Carlo (iHMC) has large errors for higher values of m, but our method performs consistently across a wide range of massimbalances. We also show unpublished data from the CL approach, shared with us by the authors of Ref. [21]. These CL results seem to converge to the correct values, in contrast to the repulsive case discussed earlier. The  10. The two plots shown above compare the systematic errors in extracting the ground state energy E 0 N ↑ ,N ↓ due to the fitting range in β used to fit the data to Eq. (33). The data is from the mass-balanced case (m = 0) with N ↑ = N ↓ = 3 particles with γ = −3.0 and LX = 40 (see Fig. 9). The left plot shows the fit for the range 25 ≤ β ≤ 75, while the plot on the right is for the fit in the range 100 ≤ β ≤ 500. The extracted ground state energy along with the error is shown as a shaded region. The exact answer from the Bethe ansatz is shown at a dotted line in the both the plots.
prediction from second order perturbation theory is also plotted. At m = 0, the exact Bethe ansatz result clearly disagrees with our algorithm. However, our results are obtained at β = 100 and the exact answer is entirely within the O(1/β) error expected from higher excited states. We discuss this issue in more detail in the next section.

VII. EXTRACTING THE GROUND STATE ENERGY
In the above section, we presented results for E at a fixed value of β = 100. We found in Fig. 6 that the results almost agreed with the exact ground state. The small disagreement was shown to be an effect of β not being large enough in Fig. 7. This suggests that β = 100 is not guaranteed to be large in all the cases. In fact, in Fig. 9 we found that our value for the energy obtained at β = 100 disagreed with the exact ground state energy from the Bethe ansatz at m = 0 by several standard deviations. This shows that it is important to be able to perform a systematic extrapolation to the β → ∞ limit to extract the ground state energy.
At m = 0.0, Fig. 9 shows a deviation of roughly 0.009 between the exact energy and the energy at β = 100 in bare units (with E FG = 0.056). While this is within the expected O(1/β) corrections, it is important to be able to perform a systematic extrapolation to β → ∞ to extract the ground-state energy. The traditional procedure followed in the literature (which we refer to as Method I) is to compute E at several values of β and then fit to the form While this approach is surely reasonable at sufficiently large values of β, it is a priori not clear what range of β should be chosen for the fit. We believe a common misconception is that the correct range of β can be determined by increasing the upper limit of β until the above form begins to fit the data well in a range. This of course depends on the precision to which the average energies are computed. Here we show that even if the errors are in the one-percent range, which is usually difficult in many cases, we can get a good fit but with wrong results if we do not choose a sufficiently large range of β.
To demonstrate the problem, in Fig. 10 we show the fit for two ranges of β at γ = −3.0 in the mass-symmetric case where we know the exact answer, shown as the dashed line. As can be seen from the figure, the fit in the range β ∈ [40, 80] (β F ∈ [1.11, 2.22]) is excellent but gives us E 0 N ↑ ,N ↓ /E FG = −2.285 (19), which is different from the exact result of E 0 N ↑ ,N ↓ /E FG = −2.40836 by several standard deviations. On the other hand notice that the fit in the range β ∈ [70, 500] (β F ∈ [1.94, 13.88]) gives a better answer of E 0 N ↑ ,N ↓ /E FG = −2.4104 (30). This also explains the deviation in Fig. 9 at a fixed value of β = 100 noted in the previous paragraph. way to compute E 0 N ↑ ,N ↓ , which can be compared with the above method. Below we discuss one such method (which we refer to as Method II) which we believe gives better precision although the analysis is more involved. Since the worm algorithm allows us to study a variety of particle number sectors efficiently, we can efficiently build the desired particle-number sector by adding one particle at a time. Consider a situation where we can tune the chemical potentials µ σ close to a critical value so that the average particle numbers fluctuate between (N ↑ , N ↓ ) and (N ↑ +1, N ↓ +1). Near such a critical point, we can develop efficient worm algorithms to sample configurations where N ↑ fluctuates by one while N ↓ remains fixed. We can then accurately measure the ratios like where Z N ↑ ,N ↓ µ as the partition function defined in Eq. (5) restricted to a fixed particle number sector, Assuming β is sufficiently large so that only the ground states contribute, we must have where is the difference in the ground-state energies of the two particle number sectors, and g is the ratio of their degeneracies, which is a fraction typically made up of small integers. If g can be determined from the knowledge that it is a fraction of small integers, the fit function (36) has a single free parameter (µ c ) for all values of β (sufficiently large) and µ ↑ . This fitting procedure gives very precise values for the critical chemical potential µ c . A similar procedure can be adapted to compute the dif- Absolute energies can be computed by adding such differences.
In order to demonstrate that method II gives us more accurate answers as compared to method I and more importantly does not depend on the range of β used in the analysis we have extracted the ground state using both the methods at m = 0.5 and γ = −0.3. However, for this study we limited our effort to a fixed non-zero value of ε = 0.01. It should be noted that when ε = 0 the definition of the ground state energy obtained from the two methods can disagree due to O(ε) errors. But it can still give us a sense of the magnitude of the errors in computation. The extrapolation to ε → 0 requires further effort but can be done if necessary.
In Fig. 11 we show our results for method II. We begin with the system containing a single particle (N ↑ = 1, N ↓ = 0) whose ground state energy is known to be zero. We then add particles slowly to reach N ↑ = N ↓ = 3 in five steps. The results from each step are shown in the figure by plotting the ratio R σ N ↑ ,N ↓ as a function of µ. The solid lines in each plot are a combined fit to the form Eq. (36) with one parameter µ c . This gives us the following results: (35).
Adding the all the values of µ c we obtain E 0 3,3 = −0.089016(53) which gives E 0 3,3 /E FG = −1.60341(95). To compare the errors obtained from method II with method I we compute E at several values of β. In Fig. 12 we show fits of this data to the form Eq. (33) in two different ranges of β. Again we see that in method I, the range of β at larger values gives a slightly different estimate of the ground-state energy as compared to the range at smaller values. Due to time discretization errors that are present at ε = 0.01, the estimate of the ground state energy computed using the two methods can disagree. Here we focus on the error in this estimate and observe that it is a factor of three more in method I as compared to method II, although individual data points were all obtained with the same precision of roughly one percent.

VIII. HIGHER DIMENSIONS
The quantum Monte Carlo method we have developed in this work is guaranteed to work well only for systems with hard-core bosons and not fermions. In one dimension, where fermions are identical to hard-core bosons in many cases, our method can be used to study fermionic systems as well. We have demonstrated this in Section VI. However, as we pointed out in Section V, we can also study a system of fermions in higher dimensions if we can compute the fermion sign S accurately. We believe this may be possible in the context of few-body physics by combining ideas of fermion bags proposed recently [50] with ideas of exponential error reduction proposed several years ago [51]. While a complete study of this more difficult problem is a topic for another paper, here we provide some evidence that our method can indeed be used for computing the ground-state energy even in higher dimensions with fermions. For this purpose, we compute the ground state energy of our Hamiltonian (1) in three spatial dimensions with L X = 8 lattice for the mass-balanced system (m ↑ = m ↓ = 1) with N ↑ = N ↓ = 2 at U = −3.9570, which corresponds to the unitary fixed point in the continuum limit. The exact ground-state energy for this system was computed some years ago as a benchmark calculation by an explicit diagonalization of the Hamiltonian using the Lanczos algorithm and was found to be E 0f 2,2 ≈ 0.1042 [44]. Here we reproduce this result using our approach. We first performed calculations at several values of β ranging from 10 to 50 at ε = 0.01. The results for the bosonic energy E are shown in Table III. Since the bosonic energy does not change much between β = 30 and 50 we assume that β = 30 is sufficiently large and perform a careful extrapolation of the ε errors there. In Fig. 13 we show our results for this extrapolation for both bosonic energy E and the fermion sign S . Note that errors due to a non-zero value of ε are linear at leading order and hence can be large. The solid lines in the figure are fits that we used to extract ε → 0 limit. In this time continuum limit we find that E = −0.1675 (20) and S = 0.00025 (2). Using Eq. (24) we estimate E 0f 2,2 ≈ 0.1090 (40), which is in reasonable agreement with the exact result E 0f 2,2 = 0.1042 [44].

IX. CONCLUSIONS
In this work, we proposed a worldline based approach to few-body physics where fermions are formulated as hard-core bosons, since they incorporate one of the ingredients of the Pauli exclusion principle, which is that two identical fermions cannot exist at the same spacetime point. On the other hand, hard-core bosons do not capture the fermion permutation sign, which needs to be taken into account explicitly before our method can truly be applicable to fermionic systems. Fortunately, in one spatial dimension, fermion permutations can only occur over the boundaries and the fermion permutation sign is positive in many cases.
Our approach is complementary to the well-developed AFQMC methods for fermionic systems, which unfortunately suffer from sign problems even in one spatial dimension. To demonstrate the power of our method, we showed that we can reproduce some exact results for the one-dimensional Hubbard model, obtained using the Bethe ansatz, for a wide range of couplings in the massbalanced case. Our method can easily be applied to the mass-imbalanced case where a general exact solution is not known. We used our approach to show that the results from the CL method, recently presented in Ref. [21], yields wrong values for repulsive interactions. This must be related to the 'fat-tailed' distributions of the observables in the CL method, as noted in the appendix of Ref. [21]. On the other hand, based on the data shared with us by the authors of that paper, the CL method seems more robust on the attractive side in the parameter range studied.
Extending our approach to higher dimensions is straightforward, although we have to confront the fermion sign problem which is equally severe for all types of interactions. Unlike the AFQMC methods, there is no particular advantage for attractive interactions as compared to repulsive interactions, although we do see that the sign problem becomes slightly milder in the attractive case. We presented some evidence that, at least for few-body physics, we may be able tame the sign problem using ideas of fermion bags and exponential error reduction. In particular, we were able to reproduce a benchmark calculation done a few years ago with four fermions at unitarity [44].
Finally, unlike the AFQMC method, our worldline approach can be formulated directly in the time continuum limit [45,46,50], although in this work we did not exploit this advantage. Instead, we studied the time discretization errors. We found that they are linear in the temporal lattice spacing ε and can be eliminated by a simple extrapolation. Without such an extrapolation we would not have been able to compare our results with the exact results from the Bethe ansatz.

Discrete time results
We first verify that our algorithm is able to reproduce the exact partition function on a small space-time lattice where we can enumerate all configurations. This should also help clarify the definition of the finite-ε transfer matrix that we use. The exact expression for the finite-ε partition function is given by (see Eq. (12)) where C are the worldline configurations. On a 2×2 spacetime lattice we can explicitly enumerate all configurations C, which gives us where Z σ 0 = 1 + 4(W ↑ m ) 2 + 4(W ↑ m ) 4 , and the local weights W σ m , W σ h , W I were defined in Section III to be Recall that the weight for each configuration is given by Eq. (17), so the number of particles N σ , number of hops n σ h and number of interactions n I for each configuration can simply be read off from the exponents in expression (B1) above. The average energy is then computed using the definition (19). Table IV compares the results for the average particle numbers and average energy obtained from our Monte Carlo method with those from the exact expression in Eq. (B1).

Continuous time results
Since our time discretization procedure is very different from conventional approaches, it is useful to verify that our results do agree with the partition function (5) in the ε → 0 limit. We can do this for small spatial lattices by explicitly diagonalizing the Hamiltonian. The dimension of the relevant Hilbert space for (N ↑ , N ↓ ) hard-core bosons in a d-dimensional spatial lattice of size L X is L d To compare with the exact Hamiltonian results, we take the ε → 0 limit by performing calculations at several values of ε and performing a linear extrapolation. Table V shows a comparison of our Monte Carlo data with exact results at a few sets of parameters in d = 1, 2, 3 dimensions for = 0.001, 0.0005 and → 0. We find that To demonstrate that extrapolation to the ε → 0 limit is necessary for the agreement, we show an illustrative fit for the fourth row in Figure 14.
a proper extrapolation in higher dimensions is important to reproduce the exact results within errors. We show an example of this in Fig. 14. We believe this agreement provides another non-trivial check of our approach in arbitrary dimensions, coupling strengths, mass-imbalances and particle numbers.

Appendix B: Second-Order Perturbation Theory
In many figures, we show results from first-and secondorder perturbation theory. Here we provide the expressions used to obtain these results. We assume the particles to be fermions for the Hamiltonian in Eq. (1), and use  Table V.
fermions in a box of size L X (even) with periodic boundary conditions.
The energy spectrum for the free theory (U = 0) can be constructed using single particle energy eigenvalues which are characterized by integers k σ = (k σ,1 , . . . , k σ,Nσ ) for σ = ↑, ↓, which label the energy eigenstates of individual fermions of each type. Let the single σ-particle energy corresponding to the integer k be t σ (k) where (k) = 4 sin 2 πk L X k = 0, ±1, . . . , ±L X /2. (B1) The total energy of a system of (N ↑ , N ↓ ) non-interacting fermions is then When N σ is odd (the values which we consider in this paper) the ground state is unique and corresponds to the choice k 0 σ ≡ (−(N σ − 1)/2, −(N σ − 1)/2 + 1, . . . , (N σ − 1)/2). The ground state energy of the full Hamiltonian in perturbation theory up to second order is then given by where the primed sum is over states (k ↑ , k ↓ ) that are related to the free ground state by an excitation of exactly one particle of each type, such that the total change in momentum is zero: k ↑ + k ↓ = k 0 ↑ + k 0 ↓ .