\section{Generalized Metropolis with force bias \label{appendix:fb}}

In this appendix, we describe our second approach using  the
generalized 
Metropolis procedure to accelerate the sampling of paths in auxiliary field (AF) space. 
We introduce a dynamic force bias, analogous to what is employed in the branching random walk methods
in constrained path or phase-free AFQMC \cite{AFQMC-lecture-notes-2013}, in proposing the updates of the 
field values, which improves the acceptance ratio and hence the MC efficiency.

 To facilitate the description of the sampling algorithm
 we first give a brief sketch of the standard path-integral AFQMC approach, on which more detailed descriptions
can be found in, for example, Refs.~\cite{assaad-lec-note} and \cite{AFQMC-lecture-notes-2013}.
Ground state AFQMC measures the static properties by
\begin{equation}
\label{eqn:projection}
 \langle \hat{O}\rangle =\frac{\langle \psi_T|\,\exp(-\beta \hat{H}/2)\,\hat{O}\,\exp(-\beta\hat{H}/2)\,|\psi_T \rangle}{ \langle \psi_T|\,\exp(-\beta \hat{H})\,|\psi_T \rangle}\,,
\end{equation}
where the Hamiltonian $\hat{H}\equiv \hat{K}+\hat{V}$ is given by Eq.~(\ref{eq:H}). 
We apply the usual Trotter-Suzuki breakup 
\begin{equation}
 e^{-\Delta \tau \hat{H}} \backsimeq e^{-\Delta \tau \hat{K}/2} e^{-\Delta \tau \hat{V}} e^{-\Delta \tau \hat{K}/2}\,
\end{equation}
and the Hubbard-Stratonovich (HS) decomposition 
\cite{PhysRevB.28.4059}
\begin{align}
 e^{\Delta \tau U n_{i\uparrow}n_{i\downarrow}} 
 &=\frac{1}{2} \sum_{x_i=\pm1} e^{(\gamma x_i -\Delta \tau U/2) (n_{i\uparrow} + n_{i\downarrow}-1)} 
 \label{eq:HS-Hirsch-charge} \\
 &\equiv 
\frac{1}{2} \sum_{x_i=\pm1} \hat{b}_i(x_i)\,, \nonumber  
\end{align}
with $ \cosh(\gamma) = \exp( -\Delta \tau U/2 )$,  
arriving at the form 
\begin{equation}
      e^{-\Delta\tau {\hat H}}=
\int d{\mathbf x}\,  p({\mathbf x}) {\hat B}({\mathbf x})\,,
\label{eq:HS}
\end{equation}
where ${\mathbf x}=\{x_1, x_2, \cdots,x_{{\mathcal N}_s}\}$
The probability density function 
 $p({\mathbf x})$ is uniform for the $2^{{\mathcal N}_s}$ AF configurations under the choice of HS in Eq.~(\ref{eq:HS-Hirsch-charge}), and the one-body propagator
 is ${\hat B}({\mathbf x})= e^{-\Delta \tau \hat{K}/2} \prod_i 
 \hat{b}_i(x_i)\,e^{-\Delta \tau \hat{K}/2} $.
 
  The expression in Eq.~(\ref{eqn:projection}) is then re-written as a path integral of $M\equiv \beta/\Delta\tau$ 
 time slices. Let us consider the $l$-th time slice, and introduce the notation
 \begin{eqnarray*}
\langle\psi_l| & = &
\langle \psi_T|\,
{\hat B}({\mathbf x}^{(M)}){\hat B}({\mathbf x}^{(M-1)}) \cdots
{\hat B}({\mathbf x}^{(l+1)})\,e^{-\Delta \tau \hat{K}/2} \\
|\psi_r\rangle & = &
e^{-\Delta \tau \hat{K}/2}\,{\hat B}({\mathbf x}^{(l-1)}){\hat B}({\mathbf x}^{(l-2)}) \cdots {\hat B}({\mathbf x}^{(1)})
\,|\psi_T\rangle\,,
\end{eqnarray*}
which are both single Slater determinant wave functions if we choose $|\psi_T\rangle$ to be a Slater determinant. 
The integrand of the path-integral in the denominator of Eq.~(\ref{eqn:projection}) then becomes 
\begin{equation}
{\mathcal W}({\mathbf x})=p({\mathbf x})\, \langle \psi_l |  \prod_{i=1}^{{\mathcal N}_s}  \hat{b}_i(x_i) |\psi_r\rangle\,,
\end{equation}
where ${\mathbf x}$ denotes the collection of AF at time slice $l$.
In the standard way of sampling ${\mathcal W}$, one 
proposes to flip each auxiliary-field $x_i$ one by one, and sweeps through ${\mathbf x}$. We will update 
the entire configuration ${\mathbf x}$ (or a sub-cluster of ${\mathbf x}$ for very large system sizes), simultaneously. We define a force bias  \cite{AFQMC-lecture-notes-2013}:
\begin{equation}
 \bar{n}_{i\sigma}= \frac{ \langle \psi_l |n_{i\sigma} |\psi_r\rangle  } {\langle \psi_l  |\psi_r\rangle}\,,
 \label{eq:FB-dynamic}
\end{equation}
and propose updates of the fields with the probability density: 
\begin{equation}
{\mathcal P}({\mathbf x})\propto p({\mathbf x}) \prod_{i=1}^{{\mathcal N}_s} e^{\gamma x_i  (\bar{n}_{i\uparrow} + \bar{n}_{i\downarrow}-1)}\, 
\end{equation}
which can be sampled directly. Detailed balance then leads to a Metropolis acceptance probability given by
\begin{equation}
  {\mathcal A}({\mathbf x}\rightarrow {\mathbf x}') = \min\{1,\frac{{\mathcal W}({\mathbf x}^{\prime})\,{\mathcal P}({\mathbf x})}{{\mathcal W}({\mathbf x})\,{\mathcal P}({\mathbf x}^{\prime})}\}\,.
\end{equation}
Note that the probability function for proposing transitions
does not depend on the ``current'' configuration of AF, i.e., ${\mathcal P}({\mathbf x}\rightarrow {\mathbf x}')={\mathcal P}({\mathbf x}')$.
If ${\mathcal P = W}$, all updates will be accepted. 
Because of the force bias, ${\mathcal P}$ approximates ${\mathcal W}$ up to ${\mathcal O}(\sqrt{\Delta \tau})$, 
leading to typically high acceptance ratio.

Although we have used the discrete charge HS decomposition, 
the algorithm generalizes straightforwardly 
to continuous HS transformations. 
We comment that the use of the dynamic force bias in Eq.~(\ref{eq:FB-dynamic}) effectively introduces a 
background subtraction \cite{PhysRevB.88.125132,AFQMC-lecture-notes-2013} in the decomposition of  Eq.~(\ref{eq:HS-Hirsch-charge}). 
That is, if one were to employ the standard updating algorithms
\emph{without the force bias}, one would find Eq.~(\ref{eq:HS-Hirsch-charge})
much less efficient than a continuous charge decomposition which subtracts a constant 
background. This discrepancy in efficiency grows more as the system density decreases, which is especially
relevant since
the systems studied here are at the low density limit. (See Ref.~\cite{PhysRevB.88.125132}
for an analysis of the efficiency of HS transformations, and Ref.~\cite{PhysRevE.70.056702} for 
discussion on how the dynamic force bias automatically introduces an optimal constant background shift.)

Some other features of our algorithm are:

\begin{itemize} 

\item Since we always work in the dilute limit, the memory is saved by only storing the wave function and calculating the Green function on the fly. We divide the path of $M$ slices into 
$\sqrt{M}$ blocks, and
only track one block each time. The wave function at the beginning of each block is stored. The largest number of wave functions stored in our code is $\sim2\sqrt{M}$.

\item The wave function is transformed
between real and momentum space by fast Fourier transformation, so that all the one-body operators during projection are diagonal, and Green functions in different space are easily obtained.

\item When we only need the energy, we separate it into kinetic and potential energy.  They are diagonal either in momentum or real space,  where we do not need to calculated the whole Green function.
To improve statistics, we measure the energy
anywhere along the path  and combine them, including the mixed estimator on both side.

\item The standard determinantal QMC formalism as sketched above turns out to have a divergence of the Monte Carlo variance. We discuss the variance problem and its solution 
separately elsewhere \cite{HShi-SZhang-to-be-published}. The solution involves the introduction of a bridge link, which we 
have implemented in the calculations presented here.
The force bias and basic sampling algorithm described above
remain unchanged.


\end{itemize}


\section{Extrapolation to the continuum limit \label{appendix:fss}}
We have described the extrapolation procedure of our lattice results to the continuum limit, and the subsequent 
analysis to reach the thermodynamic limit. Here
we illustrate the finite  size extrapolation in few-body systems.

The extrapolation to the continuum limit, 
for a fixed number of particles, 
must be consistent and independent of 
the type of kinetic energy dispersion.
For a two-body problem on the lattice, exact results can be obtained for large system sizes by 
mapping to a one-body problem in the center of mass system.
The results are shown in Fig.~\ref{fig:fss}(a), which 
fit well a 4th-order polynomial function in $1/L$. We see from the inset that the coefficient on the linear term  is zero within numerical precision.

We also show the finite size effect in the four-body problem from QMC, in Fig.~\ref{fig:fss}(b), reaching large lattice sizes. The same general behavior is seen as in the two-body problem.  
We have also studied the finite-size behavior of the BCS solution, finding similar trends but with different slopes.
In the many-body system, our QMC data are consistent with these observations as well. They are thus fitted 
 with a 4th-order polynomial function with a vanishing $1/L$ coefficient, as described in the main text.

\begin{figure*}[htbp]
  \includegraphics[scale=0.5]{Size-scaling.eps}
  \caption{\label{fig:fss} (Color online) 
  Extrapolation of finite-size lattice to the continuum limit in few body problems. Panel (a) shows 
  exact diagonalization results for the two-body problem at $\ln(a k_F)=0.5$, while panel (b) shows
  QMC solutions for the four-body problem at $\ln(a k_F)=0.0$.
In each case, results are obtained for both the Hubbard and the quadratic dispersions.
A 4th-order polynomial function in $1/L$ fits well both dispersions, and the extrapolated results in continuum limit agree well with each other. The insets indicate that the coefficients on $1/L$ are negligible in both cases.
  }
\end{figure*}



\section{Equation of state data \label{appendix:eos}}
We list the data for the equation of state in Fig.~\ref{fig:energy}.  The QMC energy data are calculated by our branching random walk approach with BCS trial wave functions.

\begin{table}
\caption{Data of the equation of state in Fig.~\ref{fig:energy}. 
The interaction strength, given in the first column, are $\ln(a k_F)=y+\ln(2)/2$, with $y$ from 
$-0.75$ to $6$ (in increments of $0.25$ up to $y=2$, then increments of 0.5 up to $y=5$).
}. 
\begin{ruledtabular}
\begin{tabular}{cccc}
$\ln(a k_F)$ & $E_{QMC}/E_{FG}$ & Error bar & $E_{BCS}/E_{FG}$\\
\hline
   -0.403426     &    -5.512634      &      0.000619    &    -4.651252  \\
   -0.153426     &    -3.262997      &      0.000487    &    -2.427641  \\
    0.096574     &    -1.884889      &      0.000325    &    -1.078969  \\
    0.346574     &    -1.027841      &      0.000453    &    -0.260958  \\
    0.596574     &    -0.487667      &      0.000335    &     0.235190  \\
    0.846574     &    -0.137058      &      0.000272    &     0.536119  \\
    1.096574     &     0.096228      &      0.000203    &     0.718642  \\
    1.346574     &     0.256943      &      0.000167    &     0.829348  \\
    1.596574     &     0.371799      &      0.000162    &     0.896491  \\
    1.846574     &     0.456471      &      0.000141    &     0.937204  \\
    2.096574     &     0.521804      &      0.000173    &     0.961859  \\
    2.346574     &     0.572904      &      0.000111    &     0.976740  \\
    2.846574     &     0.647340      &      0.000103    &     0.990927  \\
    3.346574     &     0.700067      &      0.000067    &     0.997096  \\
    3.846574     &     0.737144      &      0.000128    &     0.997307  \\
    4.346574     &     0.767283      &      0.000099    &     0.997765  \\
    4.846574     &     0.793547      &      0.000068    &     1.000206  \\
    5.346574     &     0.813073      &      0.000053    &     1.000436  \\
    6.346574     &     0.842689      &      0.000036    &     1.000654 
\end{tabular}
\end{ruledtabular}
\end{table}

\section{Full digits for Table \ref{tab:fitparam} \label{appendix:fitparam}}
\begin{table}
\caption{\label{tab:fitparamfull} 
Final parameter values (full digits) in the parametrization [Eqs.~(\ref{eqn:fitcross}-\ref{eqn:fitbec})] of the exact EOS from QMC. 
}
\begin{ruledtabular}
\begin{tabular}{cccc}
&  $a^l_0$ &  -11.804127317953723 & \\
&  $a^l_1$ &  14.675499370762239  &\\
&  $a^l_2$ &  -4.855080880566919  &\\
&  $a_0$   &  -0.819842357425408  &\\
&  $a_1$   &  0.1273251139440354  &\\
&  $a_2$   &  0.0685123559420463  &\\
&  $a_3$   &  -0.014505432043856327&\\
&  $a_4$   &  -0.009191602440101383 &\\
&  $a_5$   &  0.004190575139056055   &\\
&  $a_6$   &  -0.0006367374265820822   &\\
&  $a_7$   &  3.431232818866204$\times 10^{-5}$  &\\
&  $a^r_2$ & -0.060852400057644876  &\\
&  $a^r_3$ &  0.36401186693517423   &\\
&  $a^r_4$ & -0.61531422724189
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{acknowledgments}
We thank J.~Carlson 
for useful discussions. 
This research was supported by DOE (grant no.~DE-SC0008627), NSF 
(grant no.~DMR-1409510), and the Simons Foundation. Computing was carried out at the 
Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory,
 which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725,
and at the computational facilities at the College of William and Mary.
\end{acknowledgments}

