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

Hersh Singh; Shailesh Chandrasekharan (Duke University)

arXiv:1812.05080·hep-lat·May 1, 2019

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

Hersh Singh, Shailesh Chandrasekharan (Duke University)

PDF

TL;DR

This paper develops a worldline lattice approach with worm algorithms for simulating few-body non-relativistic bosons and fermions, revealing advantages over traditional methods and potential ways to mitigate the fermion sign problem.

Contribution

It introduces a novel worldline lattice method with worm algorithms for non-relativistic particles, enabling efficient sampling and fermion sign analysis, and assesses its limitations and potential in higher dimensions.

Findings

01

Efficient sampling of worldline configurations in fixed particle-number sectors.

02

Fermion permutation sign as an observable for energy extraction.

03

Demonstrates the method's advantages over auxiliary field approaches in 1D.

Abstract

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…

Tables5

Table 1. Table 1: Table of transition probabilities among the group of three configurations in Fig. 4 b

Case I: $W_{2} > W_{1} + W_{3}$
$P_{21} = W_{1} / W_{2}$ , $P_{23} = W_{3} / W_{2}$ , $P_{12} = P_{32} = 1$
Case II: $W_{1} > W_{2} + W_{3}$
$P_{12} = W_{2} / W_{1}$ , $P_{13} = W_{3} / W_{1}$ , $P_{21} = P_{31} = 1$
Case III: $W_{a} \leq W_{b} + W_{c}$
$P_{12} = (W_{1} + W_{2} - W_{3}) / 2 W_{1}$	$P_{13} = (W_{1} + W_{3} - W_{2}) / 2 W_{1}$
$P_{21} = (W_{2} + W_{1} - W_{3}) / 2 W_{2}$	$P_{23} = (W_{2} + W_{3} - W_{1}) / 2 W_{2}$
$P_{31} = (W_{3} + W_{1} - W_{2}) / 2 W_{3}$	$P_{32} = (W_{3} + W_{2} - W_{1}) / 2 W_{3}$

Table 2. Table 2: Table of transition probabilities among the group of four configurations in Figs. 3 b and 4 c

Case I: $W_{2} > W_{1} + W_{3} + W_{4}$
$P_{21} = W_{1} / W_{2}$ , $P_{23} = W_{3} / W_{2}$ , $P_{24} = W_{4} / W_{2}$
$P_{12} = P_{32} = P_{42} = 1$
Case II: $W_{1} > W_{2} + W_{3} + W_{4}$
$P_{12} = W_{2} / W_{1}$ , $P_{13} = W_{3} / W_{1}$ , $P_{14} = W_{4} / W_{1}$
$P_{21} = P_{31} = P_{41} = 1$
Case III: $W_{a} \leq W_{b} + W_{c}$
$P_{12} = (W_{1} + W_{2} - 2 W_{h}) / 2 W_{1}$	$P_{13} = (W_{1} + 2 W_{h} - W_{2}) / 4 W_{1}$
$P_{14} = (W_{1} + 2 W_{h} - W_{2}) / 4 W_{1}$	$P_{21} = (W_{2} + W_{1} - 2 W_{h}) / 2 W_{2}$
$P_{23} = (W_{2} + 2 W_{h} - W_{1}) / 4 W_{2}$	$P_{24} = (W_{2} + 2 W_{h} - W_{1}) / 4 W_{3}$
$P_{31} = (2 W_{h} + W_{1} - W_{2}) / 4 W_{h}$	$P_{32} = (2 W_{h} + W_{2} - W_{1}) / 4 W_{h}$
$P_{41} = (2 W_{h} + W_{1} - W_{2}) / 4 W_{h}$	$P_{42} = (2 W_{h} + W_{2} - W_{1}) / 4 W_{h}$
$P_{34} = P_{43} = 1 / 2$

Table 3. Table 3: ⟨ E ⟩ delimited-⟨⟩ 𝐸 \langle E\rangle for various values of β 𝛽 \beta in the three dimensional mass-balanced model ( m ↑ = m ↓ = 1 subscript 𝑚 ↑ subscript 𝑚 ↓ 1 m_{\uparrow}=m_{\downarrow}=1 ) at L X = 8 subscript 𝐿 𝑋 8 L_{X}=8 , U = − 3.9570 𝑈 3.9570 U=-3.9570 and ε = 0.01 𝜀 0.01 \varepsilon=0.01 for N ↑ = N ↓ = 2 subscript 𝑁 ↑ subscript 𝑁 ↓ 2 N_{\uparrow}=N_{\downarrow}=2 particles.

$β$	$10$	$15$	$20$	$30$	$50$
$⟨ E ⟩$	$- 0.12 (1)$	$- 0.322 (5)$	$- 0.368 (5)$	$- 0.387 (2)$	$- 0.387 (2)$

Table 4. Table 4: Comparison of our Monte Carlo method with the exact results for a 2 × 2 2 2 2\times 2 space-time lattice in 1 + 1 1 1 1+1 dimensions. We show results for the the average energy ⟨ E ⟩ delimited-⟨⟩ 𝐸 \langle E\rangle and average particle numbers N ↑ , N ↓ subscript 𝑁 ↑ subscript 𝑁 ↓ N_{\uparrow},N_{\downarrow} .

$μ_{↑}$	$μ_{↓}$	$m_{↑}$	$m_{↓}$	$U$	$⟨ E ⟩$ $\begin{matrix} MC \\ Exact \end{matrix}$	$⟨ N_{↑} ⟩$ $\begin{matrix} MC \\ Exact \end{matrix}$	$⟨ N_{↓} ⟩$ $\begin{matrix} MC \\ Exact \end{matrix}$
$1.00$	$- 2.00$	$1.50$	$0.50$	$- 4.00$	$- 0.552 802 (53)$	$1.415 660 (25)$	$0.572 415 (24)$
					$- 0.552 791$	$1.415 684$	$0.572 407$
$2.00$	$0.00$	$1.50$	$0.50$	$- 2.00$	$0.674 356 (26)$	$1.706 027 (29)$	$0.732 806 (31)$
					$0.674 351$	$1.706 015$	$0.732 805$
$1.00$	$4.00$	$1.50$	$0.50$	$0.00$	$3.525 805 (51)$	$1.135 221 (19)$	$1.605 845 (16)$
					$3.525 821$	$1.135 242$	$1.605 851$
$2.50$	$2.00$	$1.50$	$0.50$	$2.00$	$1.739 032 (49)$	$1.574 851 (25)$	$0.333 418 (20)$
					$1.738 968$	$1.574 855$	$0.333 407$
$1.00$	$4.00$	$1.50$	$0.50$	$4.00$	$2.758 05 (10)$	$0.357 677 (24)$	$1.354 424 (28)$
					$2.758 091$	$0.357 684$	$1.354 429$

Table 5. Table 5: Comparison of the worm algorithm with exact results obtained by diagonalizing the Hamiltonian for small lattices in d = 1 , 2 , 3 𝑑 1 2 3 d=1,2,3 spatial dimensions. To demonstrate that extrapolation to the ε → 0 → 𝜀 0 \varepsilon\to 0 limit is necessary for the agreement, we show an illustrative fit for the fourth row in Figure 14 .

$d$	$L_{X}$	$β$	$μ_{↑}$	$μ_{↓}$	$m_{↑}$	$m_{↓}$	$U$		Worm Algorithm			Exact
$d$	$L_{X}$	$β$	$μ_{↑}$	$μ_{↓}$	$m_{↑}$	$m_{↓}$	$U$		$ε = 0.001$	$ε = 0.0005$	$ε \to 0$	Exact
								$⟨ E ⟩$	$8.294 02 (69)$	$8.290 95 (69)$	$8.2875 (17)$	$8.2881$
								$⟨ N_{↑} ⟩$	$2.845 25 (34)$	$2.845 02 (91)$	$2.845 32 (81)$	$2.8454$
$1$	$6$	$10$	$2.00$	$4.00$	$1.50$	$0.50$	$2.00$	$⟨ N_{↓} ⟩$	$3.952 96 (18)$	$3.953 07 (31)$	$3.952 64 (32)$	$3.9530$
								$⟨ E ⟩$	$3.677 66 (86)$	$3.6726 (23)$	$3.6663 (22)$	$3.6668$
								$⟨ N_{↑} ⟩$	$3.266 71 (48)$	$3.266 74 (48)$	$3.2659 (13)$	$3.2644$
$1$	$6$	$10$	$- 1.00$	$3.00$	$1.50$	$0.50$	$- 2.00$	$⟨ N_{↓} ⟩$	$5.153 96 (40)$	$5.152 93 (40)$	$5.151 34 (86)$	$5.1509$
								$⟨ E ⟩$	$0.091 51 (23)$	$0.091 79 (33)$	$0.092 99 (83)$	$0.0928$
								$⟨ N_{↑} ⟩$	$0.799 21 (42)$	$0.796 98 (24)$	$0.794 64 (19)$	$0.7953$
$2$	$2$	$10$	$0.40$	$0.30$	$1.50$	$0.50$	$2.00$	$⟨ N_{↓} ⟩$	$0.405 62 (38)$	$0.410 07 (26)$	$0.415 27 (76)$	$0.4148$
								$⟨ E ⟩$	$5.2475 (53)$	$5.2571 (37)$	$5.2661 (76)$	$5.2631$
								$⟨ N_{↑} ⟩$	$3.568 41 (94)$	$3.571 36 (67)$	$3.5741 (13)$	$3.5732$
$2$	$2$	$10$	$1.50$	$4.00$	$1.50$	$0.50$	$- 2.00$	$⟨ N_{↓} ⟩$	$2.626 65 (97)$	$2.630 11 (69)$	$2.6331 (16)$	$2.6328$
								$⟨ E ⟩$	$21.8894 (15)$	$21.9257 (15)$	$21.9610 (17)$	$21.9587$
								$⟨ N_{↑} ⟩$	$7.673 79 (40)$	$7.663 49 (41)$	$7.653 07 (75)$	$7.6524$
$3$	$2$	$10$	$4.80$	$5.00$	$1.50$	$0.50$	$2.00$	$⟨ N_{↓} ⟩$	$2.397 72 (48)$	$2.415 00 (48)$	$2.432 18 (63)$	$2.4321$
								$⟨ E ⟩$	$11.1153 (22)$	$11.1625 (21)$	$11.2111 (30)$	$11.2087$
								$⟨ N_{↑} ⟩$	$6.153 73 (26)$	$6.157 98 (26)$	$6.162 27 (40)$	$6.1619$
$3$	$2$	$10$	$2.50$	$3.50$	$1.50$	$0.50$	$- 2.00$	$⟨ N_{↓} ⟩$	$3.679 81 (43)$	$3.695 01 (84)$	$3.710 93 (45)$	$3.7104$

Equations113

H

H

+ \frac{U}{a ^{d}} i \sum N_{i, ↑} N_{i, ↓},

H_{a \to 0}

H_{a \to 0}

+ g \int d x ψ_{↑}^{†} (x) ψ_{↑} (x) ψ_{↓}^{†} (x) ψ_{↓} (x)

\displaystyle\langle E\rangle=\frac{1}{Z_{\mu}}\ \operatorname{Tr}\Big{(}\ H\ {e}^{-\beta H_{\mu}}\Big{)}

\displaystyle\langle E\rangle=\frac{1}{Z_{\mu}}\ \operatorname{Tr}\Big{(}\ H\ {e}^{-\beta H_{\mu}}\Big{)}

H_{μ} = H - i, σ \sum μ_{σ} N_{σ, i},

H_{μ} = H - i, σ \sum μ_{σ} N_{σ, i},

\displaystyle Z_{\mu}=\operatorname{Tr}\Big{(}e^{-\beta H_{\mu}}\Big{)}

\displaystyle Z_{\mu}=\operatorname{Tr}\Big{(}e^{-\beta H_{\mu}}\Big{)}

E_{N_{↑}, N_{↓}} = E_{N_{↑}, N_{↓}}^{0} - N_{↑} μ_{↑} - N_{↓} μ_{↓} .

E_{N_{↑}, N_{↓}} = E_{N_{↑}, N_{↓}}^{0} - N_{↑} μ_{↑} - N_{↓} μ_{↓} .

E_{N_{↑} + 1, N_{↓} + 1}^{0} - E_{N_{↑}, N_{↓}}^{0} = μ_{↑}^{c} + μ_{↓}^{c},

E_{N_{↑} + 1, N_{↓} + 1}^{0} - E_{N_{↑}, N_{↓}}^{0} = μ_{↑}^{c} + μ_{↓}^{c},

\displaystyle\langle N_{\sigma}\rangle=\frac{1}{Z_{\mu}}\ \operatorname{Tr}\Big{(}\sum_{i}\ N_{\sigma,i}\ e^{-\beta H_{\mu}}\Big{)},

\displaystyle\langle N_{\sigma}\rangle=\frac{1}{Z_{\mu}}\ \operatorname{Tr}\Big{(}\sum_{i}\ N_{\sigma,i}\ e^{-\beta H_{\mu}}\Big{)},

H_{d}

H_{d}

H_{h}

Z_{μ} =

Z_{μ} =

\displaystyle\operatorname{Tr}\Bigg{(}e^{-(\beta-t_{k})H_{d}}H_{h}e^{-(t_{k}-t_{k-1})H_{d}}H_{h}\cdots H_{h}e^{-t_{1}H_{d}}\Bigg{)},

Z_{μ} = C \sum Ω (C),

Z_{μ} = C \sum Ω (C),

Ω (C) = i, t \prod W^{i, t} (C),

Ω (C) = i, t \prod W^{i, t} (C),

W_{h}^{σ} = t_{σ} ε .

W_{h}^{σ} = t_{σ} ε .

\displaystyle W^{\sigma}_{m}=\exp\big{(}-\varepsilon(2dt_{\sigma}-\mu_{\sigma})\big{)}.

\displaystyle W^{\sigma}_{m}=\exp\big{(}-\varepsilon(2dt_{\sigma}-\mu_{\sigma})\big{)}.

W_{I} = exp (- ε U),

W_{I} = exp (- ε U),

Ω (C)

Ω (C)

⟨ E ⟩

⟨ E ⟩

E (C)

E (C)

⟨ N_{σ} ⟩ = \frac{1}{Z _{μ}} C \sum N_{σ} (C) Ω (C) .

⟨ N_{σ} ⟩ = \frac{1}{Z _{μ}} C \sum N_{σ} (C) Ω (C) .

P_{aa} = 1 - b \neq = a \sum P_{ab} .

P_{aa} = 1 - b \neq = a \sum P_{ab} .

Z_{μ}^{f} = C \sum S (C) Ω (C)

Z_{μ}^{f} = C \sum S (C) Ω (C)

⟨ S ⟩ = \frac{Z _{μ}^{f}}{Z _{μ}},

⟨ S ⟩ = \frac{Z _{μ}^{f}}{Z _{μ}},

⟨ E^{f} ⟩ = ⟨ E ⟩ - \frac{1}{β} ln (⟨ S ⟩),

⟨ E^{f} ⟩ = ⟨ E ⟩ - \frac{1}{β} ln (⟨ S ⟩),

E_{N_{↑}, N_{↓}}^{f 0}

E_{N_{↑}, N_{↓}}^{f 0}

m = \frac{m _{↑} + m _{↓}}{2}, \overline{m} = \frac{m _{↑} - m _{↓}}{m _{↑} + m _{↓}},

m = \frac{m _{↑} + m _{↓}}{2}, \overline{m} = \frac{m _{↑} - m _{↓}}{m _{↑} + m _{↓}},

γ = \frac{U}{n} .

γ = \frac{U}{n} .

E_{FG} = \frac{1}{3} N ϵ_{F},

E_{FG} = \frac{1}{3} N ϵ_{F},

2 π I_{i}

2 π I_{i}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

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

Hersh Singh and Shailesh Chandrasekharan

Department of Physics, Box 90305, Duke University, Durham, NC 27708, USA

Abstract

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 Lin and Meyer (2015).

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 Hammer (2018); Machleidt and Entem (2011). This chiral effective field theory ( $\chi$ -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 $\chi$ -EFT can require non-perturbative calculations. While analytic techniques based on resummations of Feynman diagrams are useful for up to three particles Bedaque et al. (1999), numerical approaches especially based on Monte Carlo methods become necessary for higher number of particles Lonardoni et al. (2018); Lähde et al. (2014).

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 Lee (2009); Drut and Nicholson (2013); Carlson et al. (2015). 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 Lähde et al. (2015). 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 Gross and Bloch (2017); Bloch (2005). Also, many interesting quantum phenomena in higher dimensions have analogues in one spatial dimension Endres (2012a); Drut et al. (2018); Nishida and Son (2010); Endres (2013); Yong and Son (2018). 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 Rammelmüller et al. (2017a, 2015); Loheac et al. (2015).

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. Rammelmüller et al. (2017b). Recently, this approach was also extended to higher dimensions Rammelmüller et al. (2018). Unfortunately, it is well known that CL methods may have uncontrolled systematic errors and can converge to wrong results Aarts et al. (2010). In fact, the authors of Ref. Rammelmüller et al. (2017b) 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 Yang (1967); Gaudin (1967); Lieb and Wu (1968); Turbiner (2016). With open boundary conditions or odd number of fermions with periodic boundary conditions, sign problems are absent in the worldline formulation Wiese (1993); Evertz (2003). Recently the two-dimensional lattice Thirring model with open boundary conditions was formulated using the worldline method Ayyar et al. (2018) to provide benchmark results to test the Lefschetz thimble approach Alexandru et al. (2017). Thus, we should also be able to test the recent CL results of Ref. Rammelmüller et al. (2017b) 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 Prokof’ev and Svistunov (2001). While worldline methods to study bosonic quantum field theories are well known by now Chandrasekharan (2008); Wolff (2010); Gattringer (2014), and have been used in several studies so far Banerjee and Chandrasekharan (2010); Bruckmann et al. (2016); de Forcrand et al. (2014); Gattringer et al. (2018a); Bruckmann et al. (2015); Gattringer et al. (2018b), their applicability to bosonic quantum many-body physics has remained relatively unexplored Batrouni et al. (2009); Endres (2012b). 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. Rammelmüller et al. (2017b). 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\sim 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. Bour et al. (2011) 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

[TABLE]

where $c^{\dagger}_{i,\sigma},c_{i,\sigma}^{\phantom{\dagger}}$ create and annihilate bosons of species, say, spin-up ( $\sigma={\uparrow}$ ) or spin-down ( $\sigma={\downarrow}$ ), at the site $i$ on a $d$ -dimensional spatial lattice and $\hat{\alpha}$ represents unit vectors in the $d$ spatial directions. The site occupation-number operators are defined as $N_{i,\sigma}=c^{\dagger}_{i,\sigma}c_{i,\sigma}^{\phantom{\dagger}}$ . The parameters $t_{\sigma}=1/(2m_{\sigma}a^{2})$ give the hopping strength of the particles in terms of the masses of the particles $m_{\sigma}$ 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\rightarrow 0$ . In $d=1$ , this limiting procedure leads to the continuum theory with a delta function interaction

[TABLE]

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 Bedaque et al. (1999). 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_{\uparrow},N_{\downarrow}}$ , which is the lowest eigenvalue of the Hamiltonian in Eq. 1 in the subspace with particle numbers $N_{\uparrow}$ and $N_{\downarrow}$ . One way to accomplish this is to simply compute the average energy

[TABLE]

at very low temperature ( $\beta\to\infty$ ), where we use the modified lattice Hamiltonian,

[TABLE]

with chemical potentials $\mu_{\sigma}$ for the two particle species, and the partition function

[TABLE]

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 $\langle E\rangle=E^{0}_{N_{\uparrow},N_{\downarrow}}$ in the limit $\beta\rightarrow\infty$ . Thus, we need a method to compute $\langle E\rangle$ reliably for large values of $\beta$ .

In the next two sections, we will find an expression for $\langle E\rangle$ 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_{\uparrow}$ and $N_{\downarrow}$ particles in the presence of a chemical potential, which is given by

[TABLE]

If we can tune $\mu_{\sigma}$ to critical values $\mu_{\sigma}^{c}$ such that $E_{N_{\uparrow}+1,N_{\downarrow}+1}=E_{N_{\uparrow},N_{\downarrow}}$ , that is

[TABLE]

then worm algorithms can sample worldline configurations in the particle-number sectors $(N_{\uparrow},N_{\downarrow})$ and $(N_{\uparrow}+1,N_{\downarrow}+1)$ very efficiently even at very large values of $\beta$ . We can monitor this by computing the average particle number,

[TABLE]

and making sure that it fluctuates between the sectors $(N_{\uparrow},N_{\downarrow})$ and $(N_{\uparrow}+1,N_{\downarrow}+1)$ even when $\beta$ is large Banerjee and Chandrasekharan (2010). 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_{\mu}=H_{d}-H_{h}$ , a sum of a diagonal term and a hopping term, where

[TABLE]

We then expand the partition function as

[TABLE]

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 Beard and Wiese (1996) and fermions Rubtsov et al. (2005). However, to develop the worm algorithm, it is convenient to discretize the time integrals by dividing $\beta$ into $L_{T}$ imaginary time steps of width $\varepsilon$ , such that $\beta=\varepsilon 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 ${\cal C}$ of both species of particles on a space-time lattice. We write this as

[TABLE]

where $\Omega({\cal C})$ is the Boltzmann weight of each worldline configuration. Figure 1 gives an illustration of ${\cal C}$ on a $1+1$ dimensional space-time lattice. In the next section, we construct a worm algorithm to update such worldline configurations. We usually perform calculations at several values of $\varepsilon$ and then extrapolate to the continuous time limit. We always find that the time-discretization errors are linear in $\varepsilon$ 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 Beard and Wiese (1996).

The Boltzmann weights $\Omega({\cal C})$ can be obtained as a product of local weights associated to each space-time lattice site,

[TABLE]

where $W^{i,t}({\cal C})$ is obtained from the local worldline configuration ${\cal C}$ at the space-time site $(i,t)$ and is a product of either $W^{\sigma}_{e},W^{\sigma}_{h}$ or $W^{\sigma}_{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^{\sigma}_{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

[TABLE]

If, instead, the particle moves forward in time, the weight is

[TABLE]

The interaction among the particle species is taken into account through the weight

[TABLE]

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

[TABLE]

where ${N_{\sigma}}$ is the number of $\sigma$ -particles in the configuration $\mathcal{C}$ , ${n_{h}^{\sigma}}$ 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(\varepsilon)$ errors, as

[TABLE]

where $E(\mathcal{C})$ is the energy of a worldline configuration ${\mathcal{C}}$ , defined as

[TABLE]

This expression for $E(\mathcal{C})$ can be derived from Eq. 17 by taking the appropriate derivative with respect to $\beta$ . The average particle number for each species defined in Eq. 8 is straightforward since worldline configurations have a well-defined particle number $N_{\sigma}(\mathcal{C})$ . Thus, we get

[TABLE]

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 ${\mathcal{C}}$ of the type shown in Fig. 1 Adams and Chandrasekharan (2003); Chandrasekharan and Jiang (2006). 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 Syljuåsen and Sandvik (2002). A bounce occurs when the local configuration does not change (i.e, the 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 three- and 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 $\varepsilon$ , 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

[TABLE]

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 $\varepsilon$ . There are three possible orderings of the weights: Case I: $W_{2}>W_{1}+W_{3}$ , Case II: $W_{1}>W_{2}+W_{3}$ and Case III: $W_{a}\leq W_{b}+W_{c}$ for all $a,b,c$ different. In Table 1, 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$ , $W_{2}=W_{m}$ (or $W_{2}=W_{m}W_{I}$ ) and $W_{3}=W_{4}=W_{h}$ . 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 1. The complete set of transition probabilities for this case is given in Table 2. We test our algorithm for small lattices in $d=1,\,2,\,3$ dimensions and the results of these tests are given in Appendix A.

V Fermions and Sign Problems

Any lattice model of hard-core bosons can be converted into a model of fermions if we take into account the fermion permutation sign. The partition function for fermions $Z^{\text{f}}_{\mu}$ , can be obtained from the bosonic one in Eq. 12 through the relation

[TABLE]

where $S(\cal C)$ is the fermion permutation sign of the worldline configuration ${\cal C}$ . If we define the average fermion permutation sign $\langle S\rangle$ as the ratio of the two partition functions in a fixed particle number sector as

[TABLE]

then the energy observable in the fermionic theory can be defined using the relation

[TABLE]

where $\langle E\rangle$ 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

[TABLE]

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, $\langle S\rangle=1$ . In the case of the open boundary conditions or in a trapping potential, the fermions cannot wind around the box and $\langle S\rangle=1$ again. This implies that the fermionic and bosonic energies are the same, that is, $E^{0f}_{N_{\uparrow},N_{\downarrow}}=E^{0}_{N_{\uparrow},N_{\downarrow}}$ . We use this result in the next two sections when we study one-dimensional problems.

In higher dimensions, we expect $\langle S\rangle<1$ . In particular, the average permutation sign will suffer from a severe signal-to-noise ratio problem when $\beta$ 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 $\langle S\rangle$ as a function of the particle number $N=N_{\uparrow}+N_{\downarrow}$ in the mass-balanced repulsive model with $m_{\sigma}=1$ and $U=4$ at temperatures such that $\beta(1-\cos(2\pi/L))\sim 1.5$ . For even $N$ we choose $N_{\uparrow}=N_{\downarrow}$ and for odd $N$ we use $N_{\uparrow}=N_{\downarrow}+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 Huffman and Chandrasekharan (2017) and exponential error-reduction techniques Luscher and Weisz (2001), 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.

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_{\uparrow}=N_{\downarrow}$ ). 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 $\langle E\rangle$ (see Eq. 3) at a fixed $\beta=100$ . In the next section, we discuss how our results change with $\beta$ 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 $\varepsilon\to 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_{\uparrow}+N_{\downarrow}$ and the particle-number density by $n=N/L_{X}$ . We define the average mass and the mass-imbalance parameter as

[TABLE]

respectively. We work in units with $m_{\sigma}=1$ and $a=1$ and parametrize the interaction strength by $\gamma$ , which is related to the bare coupling $U$ in the lattice model (1) by

[TABLE]

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_{\text{FG}}$ in the continuum defined as

[TABLE]

where $\epsilon_{\text{F}}=\pi^{2}n^{2}/8m$ is the Fermi energy.

VI.1 Mass-Balanced Systems

The continuum model (2) for fermions with $\overline{m}=0$ is known as the Gaudin-Yang model and can be solved exactly when $g>0$ using nested Bethe ansatz Yang (1967); Gaudin (1967). 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\geq 0$ , and number of particles $N=N_{\uparrow}+N_{\downarrow}$ as long as $m_{\uparrow}=m_{\downarrow}$ Lieb and Wu (1968). 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 $\overline{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_{\downarrow}=M$ spin-down fermions, and $N_{\uparrow}=N-M$ spin-up fermions. Let us normalize our energies with $t_{\uparrow}=t_{\downarrow}=t$ . Let $p_{i}=(p_{1},\dotsc,p_{N})$ and $\lambda_{\alpha}=(\lambda_{1},\dotsc,\lambda_{M})$ be two sets of ascending real numbers. Then the following $N+M$ coupled non-linear equations in the $N+M$ variables $\{p_{i},\lambda_{\alpha}\}$ determine the complete spectrum of this system:

[TABLE]

where $i=1,\dotsc,N$ , $\alpha=1,\dotsc,M$ ,

[TABLE]

and $I_{i},J_{\alpha}$ 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\sum_{i}\left(1-\cos(p_{i})\right)$ . For the ground state, we must choose $I_{i}=-\frac{1}{2}(N+1)+i$ and $J_{\alpha}=-\frac{1}{2}(M+1)+\alpha$ . The ground state for the attractive Hubbard model can be obtained from the repulsive case using the relation

[TABLE]

The exact ground-state energy obtained by this procedure for $N_{\uparrow}=N_{\downarrow}=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 $\varepsilon\to 0$ even at $\beta=100$ . For comparison, we also plot the CL results from Ref. Rammelmüller et al. (2017b) on the same plot. We notice that CL reproduces the results quite well for small values of the couplings, but starts to deviate when $\gamma\gtrsim 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 $\gamma$ 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 $\beta\epsilon_{\text{F}}=\infty$ ). We also show the worm algorithm results at different values of $\beta\epsilon_{\text{F}}$ , which indeed converge to the correct lattice answer in the limit $\beta\epsilon_{\text{F}}\to\infty$ .

VI.2 Mass-Imbalanced Systems

In contrast to the mass-balanced case, no exact solution is known for the general mass-imbalanced case ( $\overline{m}\neq 0$ ). However, our algorithm is again very efficient even in such cases. In this section, we compute some spin-balanced results and check our results against second-order perturbation theory and find good agreement for small values of $\gamma$ . Figure 8 shows our results for the ground state energy of $5+5$ particles as a function of the coupling for a high mass-imbalance of $\overline{m}=0.8$ . We perform computations at $\beta=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 $\gamma\sim 3.0$ , while the CL results disagree with both perturbation theory and our method when $\gamma\gtrsim 2$ . In the right plot of Fig. 8, we extend this to the regime of very strong repulsion. It was suggested in Ref. Rammelmüller et al. (2017b) that the flattening of the ground-state energy as a function of $\gamma$ for strong repulsive couplings, obtained using the CL method, could be a physical effect. However, the disagreement with exact Bethe ansatz calculations at $\overline{m}=0$ (Fig. 6) and with the worm algorithm at a high mass-imbalance of $\overline{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. Rammelmüller et al. (2017b) for $N_{\uparrow}=N_{\downarrow}=3$ fermions with a fixed attractive coupling $\gamma=-3.0$ across a range of mass-imbalances. We observe that imaginary-mass Hybrid Monte Carlo (iHMC) has large errors for higher values of $\overline{m}$ , but our method performs consistently across a wide range of mass-imbalances. We also show unpublished data from the CL approach, shared with us by the authors of Ref. Rammelmüller et al. (2017b). These CL results seem to converge to the correct values, in contrast to the repulsive case discussed earlier. The prediction from second order perturbation theory is also plotted. At $\overline{m}=0$ , the exact Bethe ansatz result clearly disagrees with our algorithm. However, our results are obtained at $\beta=100$ and the exact answer is entirely within the $O(1/\beta)$ 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 $\langle E\rangle$ at a fixed value of $\beta=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 $\beta$ not being large enough in Fig. 7. This suggests that $\beta=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 $\beta=100$ disagreed with the exact ground state energy from the Bethe ansatz at $\overline{m}=0$ by several standard deviations. This shows that it is important to be able to perform a systematic extrapolation to the $\beta\rightarrow\infty$ limit to extract the ground state energy.

At $\overline{m}=0.0$ , Fig. 9 shows a deviation of roughly $0.009$ between the exact energy and the energy at $\beta=100$ in bare units (with $E_{\text{FG}}=0.056$ ). While this is within the expected $O(1/\beta)$ corrections, it is important to be able to perform a systematic extrapolation to $\beta\to\infty$ to extract the ground-state energy. The traditional procedure followed in the literature (which we refer to as Method I) is to compute $\langle E\rangle$ at several values of $\beta$ and then fit to the form

[TABLE]

While this approach is surely reasonable at sufficiently large values of $\beta$ , it is a priori not clear what range of $\beta$ should be chosen for the fit. We believe a common misconception is that the correct range of $\beta$ can be determined by increasing the upper limit of $\beta$ 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 $\beta$ .

To demonstrate the problem, in Fig. 10 we show the fit for two ranges of $\beta$ at $\gamma=-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 $\beta\in[40,80]$ ( $\beta\epsilon_{\text{F}}\in[1.11,2.22]$ ) is excellent but gives us $E^{0}_{N_{\uparrow},N_{\downarrow}}/E_{\text{FG}}=-2.285(19)$ , which is different from the exact result of $E^{0}_{N_{\uparrow},N_{\downarrow}}/E_{\text{FG}}=-2.40836$ by several standard deviations. On the other hand notice that the fit in the range $\beta\in[70,500]$ ( $\beta\epsilon_{\text{F}}\in[1.94,13.88]$ ) gives a better answer of $E^{0}_{N_{\uparrow},N_{\downarrow}}/E_{\text{FG}}=-2.4104(30)$ . This also explains the deviation in Fig. 9 at a fixed value of $\beta=100$ noted in the previous paragraph.

These observations suggest a need for a complementary way to compute $E^{0}_{N_{\uparrow},N_{\downarrow}}$ , 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 $\mu_{\sigma}$ close to a critical value so that the average particle numbers fluctuate between $(N_{\uparrow},N_{\downarrow})$ and $(N_{\uparrow}+1,N_{\downarrow}+1)$ . Near such a critical point, we can develop efficient worm algorithms to sample configurations where $N_{\uparrow}$ fluctuates by one while $N_{\downarrow}$ remains fixed. We can then accurately measure the ratios like

[TABLE]

where $Z_{\mu}^{{N_{\uparrow},N_{\downarrow}}}$ as the partition function defined in Eq. 5 restricted to a fixed particle number sector,

[TABLE]

Assuming $\beta$ is sufficiently large so that only the ground states contribute, we must have

[TABLE]

where

[TABLE]

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 ( $\mu_{c}$ ) for all values of $\beta$ (sufficiently large) and $\mu_{\uparrow}$ . This fitting procedure gives very precise values for the critical chemical potential $\mu_{c}$ . A similar procedure can be adapted to compute the difference $E^{0}_{N_{\uparrow},N_{\downarrow}+1}-E^{0}_{N_{\uparrow},N_{\downarrow}}$ . 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 $\beta$ used in the analysis we have extracted the ground state using both the methods at $\overline{m}=0.5$ and $\gamma=-0.3$ . However, for this study we limited our effort to a fixed non-zero value of $\varepsilon=0.01$ . It should be noted that when $\varepsilon\neq 0$ the definition of the ground state energy obtained from the two methods can disagree due to $O(\varepsilon)$ errors. But it can still give us a sense of the magnitude of the errors in computation. The extrapolation to $\varepsilon\rightarrow 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_{\uparrow}=1,N_{\downarrow}=0$ ) whose ground state energy is known to be zero. We then add particles slowly to reach $N_{\uparrow}=N_{\downarrow}=3$ in five steps. The results from each step are shown in the figure by plotting the ratio $R^{\sigma}_{N_{\uparrow},N_{\downarrow}}$ as a function of $\mu$ . The solid lines in each plot are a combined fit to the form Eq. 36 with one parameter $\mu_{c}$ . This gives us the following results:

[TABLE]

Adding the all the values of $\mu_{c}$ we obtain $E^{0}_{3,3}=-0.089016(53)$ which gives $E^{0}_{3,3}/E_{\text{FG}}=-1.60341(95)$ .

To compare the errors obtained from method II with method I we compute $\langle E\rangle$ at several values of $\beta$ . In Fig. 12 we show fits of this data to the form Eq. 33 in two different ranges of $\beta$ . Again we see that in method I, the range of $\beta$ 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 $\varepsilon=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 $\langle S\rangle$ accurately. We believe this may be possible in the context of few-body physics by combining ideas of fermion bags proposed recently Huffman and Chandrasekharan (2017) with ideas of exponential error reduction proposed several years ago Luscher and Weisz (2001). 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_{\uparrow}=m_{\downarrow}=1$ ) with $N_{\uparrow}=N_{\downarrow}=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}\approx 0.1042$ Bour et al. (2011). Here we reproduce this result using our approach.

We first performed calculations at several values of $\beta$ ranging from $10$ to $50$ at $\varepsilon=0.01$ . The results for the bosonic energy $\langle E\rangle$ are shown in Table 3. Since the bosonic energy does not change much between $\beta=30$ and $50$ we assume that $\beta=30$ is sufficiently large and perform a careful extrapolation of the $\varepsilon$ errors there. In Fig. 13 we show our results for this extrapolation for both bosonic energy $\langle E\rangle$ and the fermion sign $\langle S\rangle$ . Note that errors due to a non-zero value of $\varepsilon$ are linear at leading order and hence can be large. The solid lines in the figure are fits that we used to extract $\varepsilon\rightarrow 0$ limit. In this time continuum limit we find that $\langle E\rangle=-0.1675(20)$ and $\langle S\rangle=0.00025(2)$ . Using Eq. 24 we estimate $E^{0f}_{2,2}\approx 0.1090(40)$ , which is in reasonable agreement with the exact result $E^{0f}_{2,2}=0.1042$ Bour et al. (2011).

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 mass-balanced 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. Rammelmüller et al. (2017b), 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. Rammelmüller et al. (2017b). 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 Bour et al. (2011).

Finally, unlike the AFQMC method, our worldline approach can be formulated directly in the time continuum limit Beard and Wiese (1996); Rubtsov et al. (2005); Huffman and Chandrasekharan (2017), 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 $\varepsilon$ 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.

Acknowledgments

We would like to thank Jens Braun, Joaquín Drut and Lukas Rammelmüller for sharing their data published in Ref. Rammelmüller et al. (2017b), the unpublished data that appears in Fig. 9, and for productive discussions about their results. We would also like to thank Dean Lee, Abhishek Mohapatra and Roxanne Springer for helpful conversations about effective field theories. The material presented here is based upon work supported by the U.S. Department of Energy, Office of Science, Nuclear Physics program under Award Numbers DE-FG02-05ER41368.

Appendix A Exact Results from Small Lattices

In this appendix we compare our Monte Carlo results with exact results on small lattices over a wide range of parameters. This study helps verify the correctness of our algorithm and provides some benchmark calculations for readers to understand our model and verify our results if necessary.

A.1 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- $\varepsilon$ transfer matrix that we use. The exact expression for the finite- $\varepsilon$ partition function is given by (see Eq. 12)

[TABLE]

where $\mathcal{C}$ are the worldline configurations. On a $2\times 2$ space-time lattice we can explicitly enumerate all configurations $\mathcal{C}$ , which gives us

[TABLE]

where $Z_{0}^{\sigma}=1+4(W_{m}^{\uparrow})^{2}+4(W_{m}^{\uparrow})^{4}$ , and the local weights $W_{m}^{\sigma}$ , $W_{h}^{\sigma}$ , $W_{I}$ were defined in Section III to be

[TABLE]

Recall that the weight for each configuration is given by Eq. 17, so the number of particles ${N_{\sigma}}$ , number of hops ${n_{h}^{\sigma}}$ and number of interactions ${n_{I}}$ for each configuration can simply be read off from the exponents in expression (A.1) above. The average energy is then computed using the definition (19). Table 4 compares the results for the average particle numbers and average energy obtained from our Monte Carlo method with those from the exact expression in Section A.1.

A.2 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 $\varepsilon\to 0$ limit. We can do this for small spatial lattices by explicitly diagonalizing the Hamiltonian. The dimension of the relevant Hilbert space for $({N_{\uparrow},N_{\downarrow}})$ hard-core bosons in a $d$ -dimensional spatial lattice of size $L_{X}$ is $\binom{L_{X}^{d}}{N_{\uparrow}}\binom{L_{X}^{d}}{N_{\downarrow}}$ where $N_{\uparrow},N_{\downarrow}\in\{0,\dotsc,L_{X}^{d}\}$ .

To compare with the exact Hamiltonian results, we take the $\varepsilon\to 0$ limit by performing calculations at several values of $\varepsilon$ and performing a linear extrapolation. Table 5 shows a comparison of our Monte Carlo data with exact results at a few sets of parameters in $d=1,2,3$ dimensions for $\epsilon=$ 0.001 $,$ 0.0005$$ and $\epsilon\to 0$ . We find that 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 second-order 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 $E_{N_{\uparrow},N_{\downarrow}}^{\text{f}\,0}$ to denote the the ground-state energy for $(N_{\uparrow},N_{\downarrow})$ 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 $\mathbf{k}_{\sigma}=(k_{\sigma,1},\dotsc,k_{\sigma,{N_{\sigma}}})$ for $\sigma={\uparrow},{\downarrow}$ , which label the energy eigenstates of individual fermions of each type. Let the single $\sigma$ -particle energy corresponding to the integer $k$ be $t_{\sigma}\epsilon(k)$ where

[TABLE]

The total energy of a system of $({N_{\uparrow}},{N_{\downarrow}})$ non-interacting fermions is then

[TABLE]

When ${N_{\sigma}}$ is odd (the values which we consider in this paper) the ground state is unique and corresponds to the choice $\mathbf{k}_{\sigma}^{0}\equiv\left(-{({N_{\sigma}}-1)/2},-{({N_{\sigma}}-1)/2}+1,\dotsc,{({N_{\sigma}}-1)/2}\right)$ . The ground state energy of the full Hamiltonian in perturbation theory up to second order is then given by

[TABLE]

where

[TABLE]

where the primed sum $\sum^{\prime}$ is over states $(\mathbf{k}_{\uparrow},\mathbf{k}_{\downarrow})$ 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: $\mathbf{k}_{\uparrow}+\mathbf{k}_{\downarrow}=\mathbf{k}^{0}_{\uparrow}+\mathbf{k}^{0}_{\downarrow}$ .

Bibliography51

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Lin and Meyer (2015) H.-W. Lin and H. Meyer, Lattice QCD for Nuclear Physics , Vol. Lecture Notes in Physics 889 (Springer, 2015).
2Hammer (2018) H. W. Hammer, Few Body Syst. 59 , 58 (2018) . · doi ↗
3Machleidt and Entem (2011) R. Machleidt and D. Entem, Physics Reports 503 , 1 (2011) . · doi ↗
4Bedaque et al. (1999) P. F. Bedaque, H.-W. Hammer, and U. van Kolck, Phys. Rev. Lett. 82 , 463 (1999) . · doi ↗
5Lonardoni et al. (2018) D. Lonardoni, J. Carlson, S. Gandolfi, J. E. Lynn, K. E. Schmidt, A. Schwenk, and X. B. Wang, Phys. Rev. Lett. 120 , 122502 (2018) . · doi ↗
6Lähde et al. (2014) T. A. Lähde, E. Epelbaum, H. Krebs, D. Lee, U.-G. Meißner, and G. Rupak, Phys. Lett. B 732 , 110 (2014) , ar Xiv:1311.0477 [nucl-th] . · doi ↗
7Lee (2009) D. Lee, Prog. Part. Nucl. Phys. 63 , 117 (2009) , ar Xiv:0804.3501 [nucl-th] . · doi ↗
8Drut and Nicholson (2013) J. E. Drut and A. N. Nicholson, J. Phys. G 40 , 043101 (2013) , ar Xiv:1208.6556 [cond-mat.stat-mech] . · doi ↗