Improved lattice operators for nonrelativistic fermions
Abstract
In this work I apply a recently proposed improvement procedure, originally conceived to reduce finite lattice spacing effects in transfer matrices for dilute Fermi systems, to tuning operators for the calculation of observables. I construct, in particular, highly improved representations for the energy and the contact, as a first step in an improvement program for finitetemperature calculations. I illustrate the effects of improvement on those quantities with a groundstate lattice calculation at unitarity.
pacs:
03.75.Ss, 05.10.Ln, 12.38.Gc, 71.10.FdI Introduction
One of the main sources of uncertainty in Monte Carlo (MC) calculations of lattice field theories stems from finite lattice spacing effects Latticebooks . This is true in particular of nonrelativistic systems at finite density, where the continuum limit is approached by taking the limit of diluteness, such that the interparticle distance (where is the Fermi momentum) is much larger than the lattice spacing . This is achieved in practice by first considering large volumes at fixed particle density, and then taking the zerodensity limit. While this process is well defined, it is cumbersome to carry out, in part because calculations in large volumes require substantial amounts of CPU time. It is therefore useful to consider alternative approaches, based on effective field theory and renormalization group concepts Lepage ; ImprovementQCD , which treat latticespacing effects by modifying the ultraviolet (UV) dynamics of the theory. In such formulations, latticespacing effects are eliminated at finite volume, even for densities that are not low by conventional standards.
In recent work, Endres et al. Endresetal ; Endresetal2 proposed a novel way to systematically reduce latticespacing effects in calculations of nonrelativistic fermions. The method enables tuning of the lattice theory to high accuracy, such that the low end of the continuum energy spectrum is reproduced for the desired values of the twobody scattering parameters in the effective range expansion:
(1) 
where is the scattering phase shift, is the scattering length, and is the effective range. The connection between the bare lattice theory (or rather its twobody spectrum) and these physical parameters is given by Lüscher’s formula Luescher , which relates the phase shift to the energy of the twobody problem in a box of side :
(2) 
where and
(3) 
where the sum is over all 3D integer vectors (the evaluation of is discussed in Appendix A), and is the Heaviside function. Throughout this work, we shall take units such that , where is the mass of the fermions.
The twobody matching condition described above completely specifies the physics of dilute systems, such as those currently realized with ultracold atomic gases (see e.g. Ref. Experiments for a review of the experimental situation). In this work we shall restrict ourselves to such systems, neglecting the additional complications that arise for instance at higher densities or for nuclear systems, where three and higherbody forces play an important role.
We shall briefly review the work of Ref. Endresetal below, but it is useful to mention the main steps underlying the method at this point. First, one writes down the transfer matrix , representing the twobody interaction via a generalized HubbardStratonovich (HS) transformation HS ; NegeleOrlandbook . The latter contains arbitrary coefficients . The matrix elements of are then computed in the subspace of two particles at vanishing total momentum. The resulting matrix is then diagonalized in the wave subspace (or rather its lattice equivalent), and the HS coefficients are tuned using an iterative algorithm such that the eigenvalues of the proposed match , where is the temporal discretization parameter and are the energies dictated by Lüscher’s formula, for the desired box size and choice of scattering parameters.
In the case of the unitary limit, where by definition
(4) 
the above procedure yields considerable improvement in the approach to the continuum. Indeed, with a single coefficient (the simplest case, discussed in Sec. II) one may tune the scattering length , but the effective range remains finite due to latticespacing artifacts. Naturally, by introducing more parameters one may tune the effectiverange expansion in Eq. (1) to higher accuracy, thus systematically eliminating the need for extrapolations to the dilute limit and leaving only finitevolume effects unaccounted for.
It is the main objective of this work to extend the approach of Ref. Endresetal to designing not only improved transfer matrices but also improved operators for the calculation of observables. This is an important step toward reducing latticespacing effects in finitetemperature calculations. While these effects were studied at unitarity in Ref. Priviteraetal for the Hubbard model using dynamical meanfield theory, careful systematic studies in fullfledged MC calculations have only recently started to appear Endresetal ; Endresetal2 ; Carlsonetal ; Leeetal and are still restricted to the ground state. Interestingly, recent studies based on functional renormalization group techniques Braun have also started to tackle the problem of extrapolating to infinite volume and particle number in the unitary limit.
In Sec. II we analyze a basic example, as a primer to reviewing the more sophisticated formalism of Ref. Endresetal , which we explain briefly in Sec. III. In Sec. IV we present the main developments of this work, with the corresponding illustrative results and conclusions appearing in Sec. V. Finally, some of the less trivial numerical issues are explained in the Appendix.
Ii A simple example
Before proceeding to a more detailed discussion, let us analyze a simple case, to fix notation as well as ideas. Consider the following lattice Hamiltonian:
(5) 
where denotes the spin projection, is the bare lattice coupling constant, and denotes the number density operator at lattice position for spin projection . The value of is tuned to the desired twobody scattering properties by solving the twobody problem on the lattice and finding the scattering amplitude, which in this case can be done analytically.
The transfer matrix is then expressed approximately in powers of the imaginary time step using a SuzukiTrotter decomposition, for instance as follows:
(6) 
where
(8)  
The kinetic energy factor in Eq. (6) is easy to treat, as it is an exponential of a onebody operator which is diagonal in momentum space. Notice that we have chosen to define the dispersion relation in momentum space, as opposed to defining it via a discrete representation of the Laplacian operator in coordinate space (which is common in Hubbardmodel type formulations).
On the other hand, the potential energy factor represents a challenge, as it is the exponential of a nontrivial twobody operator (as the original Hamiltonian). It is at this point that the HS transformation plays a central role, allowing us to exchange the complexity of the twobody operator for a path integral involving only onebody operators. Specifically, we write
where , , and is an external auxiliary field. In this way, one decouples the transfer matrix for each spin, such that we may write
(10) 
where, up to order ,
(11) 
which in the oneparticle subspace reduces to
(12) 
Notice that in Eq. (II) we have used a version of the HS transform due to Lee DeanLee , in which the auxiliary field is continuous and compact (the integral is restricted to the interval ), as opposed to more conventional versions that are discrete, or continuous but unbounded. Apart from the path integral, at this point applying the transfer matrix becomes a problem of applying a product of onebody operators, which one can easily deal with.
All of the above steps are standard in the literature of manybody MC calculations, except perhaps for the use of a “perfect” dispersion relation defined in momentum space, which has become more common only in recent years BDM ; Carlsonetal . This feature can be regarded as a basic kind of improvement, as it reduces latticespacing effects relative to Hubbardmodel approaches, where at high momenta.
This example captures the main features of most modern manyfermion MC calculations. Nevertheless, the simplicity of this approach is in some ways excessive, largely because we only have one coupling constant at our disposal. Indeed, should we desire to fix more than one coefficient in the effectiverange expansion, we would need a richer bare interaction, and a correspondingly more sophisticated HS transformation. In this simple case, to reach for instance the unitary limit, we can tune the scattering length, but we also need to consider very dilute systems to avoid the effects of finite range induced by the UV lattice cutoff , as in that case one has . Moreover, taking for example 80 particles in a lattice volume, which shall be our illustrative example below, one finds , which is not nearly as small as one would like.
As mentioned before, one can resort to alternative approaches that modify the Hamiltonian by including higherorder terms in a lowmomentum expansion, tuning the corresponding couplings to eliminate UV effects. The latter strategy is essentially what Ref. Endresetal advocates, following the spirit of the Lattice QCD program initiated by Symanzik many years ago ImprovementQCD to design improved effective actions that approach the continuum limit (see e.g. Ref. MILC_RMP ). Since it will be useful for the rest of this work, we shall briefly review the method of Ref. Endresetal in the next section.
Iii Improved transfer matrices
The work of Ref. Endresetal tackles the problem of reducing UV lattice effects by defining an improved transfer matrix, using a generalized HS transformation and a lowmomentum expansion. This is accomplished by allowing the constant to have a nontrivial momentum dependence. For simplicity, it is useful to take to be diagonal in momentum space. Following Ref. Endresetal , we expand using a set of operators as
(13) 
where one may choose to define, for convenience,
(14) 
as done in Ref. Endresetal (up to an dependent constant in front), or, as we shall use in the rest of this work,
(15) 
Notice that this last expression contains only even powers of , and is therefore an analytic function of the momentum . Both of the above choices for behave as at low momenta. For the purposes of Monte Carlo calculations, the operator can be computed once and for all at the beginning of the calculation, and is applied at each time slice using Fourier transforms.
In order to determine the coefficients , we find the explicit form of for two particles, which upon integrating the auxiliary field is given in momentum space by
(16)  
where is the lattice volume and . One cannot fail to notice that the above transfer matrix is not Galilean invariant. We elaborate on this issue in the appendices. In this work we shall take , following the steps of Ref. Endresetal . The operators in Eqs. (14) and (15) are defined to be constant above , and the kinetic energy factor is defined to vanish above that boundary, such that there is no propagation for .
Evaluating in the centerofmass frame, we have
(17) 
where and are incoming and outgoing relative momenta. By diagonalizing this expression, we may identify the eigenvalues of with , where are the twoparticle eigenvalues of the Hamiltonian we are implicitly defining.
1  0.68419  –  –  –  – 

2  0.53153  0.07896  –  –  – 
3  0.49278  0.04366  0.01807  –  – 
4  0.47217  0.03711  0.00784  0.00467  – 
5  0.45853  0.03331  0.00718  0.00132  0.00129 
One may then tune the such that this energy spectrum matches the one required by Lüscher’s formula for the lowest eigenvalues, given the desired values of the scattering parameters in Eq. (1), and using . Specifically, if we are interested in describing the unitary limit, the eigenvalues are determined by the zeros of the function in Eq. (3). The actual fitting of the can be performed iteratively, as described in detail in Ref. Endresetal . Notice, in particular, that our expression for ceases to be Hermitian when we promote to be an operator. Therefore, one needs to use rather than to diagonalize and fit the eigenvalues, as well as for the actual Monte Carlo calculations.
The results of our fits for are shown in Table 1. The quality of the improvement can be assessed by plotting as a function of , which is shown in Fig. 1 for , , and . As can be appreciated in the figure, the same effect is achieved as in Ref. Endresetal : as the order of the expansion is increased, the transfer matrix is accurately tuned to unitarity up to progressively higher momenta.
Iv Improved observables
The above procedure represents a significant step forward in mitigating latticespacing effects in MC calculations, especially considering that it requires only a small coding investment for its implementation in extant MC codes, and it results in minimal computational overhead.
A somewhat unsettling issue remains, however, particularly in connection with improving finite temperature lattice calculations, such as those of Refs. BDM ; DLT ; DLWM . Indeed, in those calculations, as well as in similar groundstate approaches, the transfer matrix is not the only object carrying latticespacing effects: the operators used to compute expectation values also suffer from the same problems. The situation may appear problematic at first sight, as the improvement strategy defined above does not directly define improved operators that we could use to compute observables at finite temperature.
On the other hand, in conventional formulations, knowledge of the explicit dependence of the transfer matrix (as in the simple example described in Sec. II) allows us to take a derivative that brings the Hamiltonian down from the exponent (c.f. Eqs. (10) and (11)), which results in a practical expression for the calculation of the energy ^{1}^{1}1 This is, of course, just a simplified case of the generating functional formalism (see e.g. Ref. ZinnJustinbook ), where one introduces a source coupled to the operator of interest, and then takes functional derivatives to bring down the desired powers of the operator from the exponent. In this example, the “source” is of course just the time step , or equivalently the inverse temperature parameter . .
To fix ideas, let us consider the grand canonical ensemble (see also Ref. BDM ; BSS ), where the partition function satisfies, by definition,
(18) 
Averages of observables can be extracted from knowledge of through derivatives, and in particular the energy is obtained by means of
(19) 
When using the formalism derived in Sec. II, we have (assuming the system is unpolarized)
(20) 
where
(21)  
(22) 
and where now is to be regarded as a spacetime varying field, and is restricted to the th imaginarytime slice. The determinant is to be taken in the subspace of oneparticle states, and
(23) 
As the derivative in Eq. (19) becomes now a derivative with respect to , it is clear that all we require is knowing how to differentiate with respect to , because
(24) 
where
(25) 
We shall take the point of view that differentiation of partition functions with respect to parameters in the Hamiltonian is the proper way to obtain expectation values of operators. We shall then generalize the improvement procedure of Ref. Endresetal to design highly improved operators (in particular for the energy and the contact) and accomplish this by taking formal derivatives of the transfer matrix and writing them in an operator expansion, fitting the expansion coefficients to reproduce the low end of the exact twoparticle spectra. By ”highly improved” we mean improvements that go well beyond tuning the first two coefficients in the effectiverange expansion.
While in the above example we have focused on the calculation of the energy, which we resume in the next section, similar derivations apply for the calculation of the contact, as we shall see in Sec. IV.2.
iv.1 Energy
In the simple case presented in Sec. II, the calculation of the first derivative results in the following expression:
(26) 
where is the following anticommutator
(27) 
and
(28) 
The derivative of with respect to appearing in Eq. (28) is known analytically in this case. When using improved transfer matrices, however, that derivative is somewhat harder to compute, as we only know the coefficients in Eq. (13) as a result of a numerical fit, which complicates the calculation of .
Here we extend the procedure of Endres et al. to fitting the coefficients in the expansion of the derivative of . There are at least two advantages in following this route over computing the derivatives via finite differences. In the first place, the fitting procedure can provide a much more accurate determination of the expansion coefficients of the derivative, with less effort. (That being said, finite differences do provide a useful starting point for the fitting routine.) Secondly, this route allows one to use different operators for different observables, regardless of what is used for itself. This feature can potentially be very convenient: different observables are in general sensitive to different regions of momentum space, such that there is no a priori reason why a given set of operators should be universally useful.
The main idea behind the method remains the same as in Ref. Endresetal . Lüscher’s formula provides us with the exact twoparticle spectrum, such that the eigenvalues of the exact twoparticle transfer matrix and its derivative(s) are known:
(29) 
where are the exact twoparticle energies in a continuous box ^{2}^{2}2 We use here for didactical purposes, but one may use the symmetric expression instead, with straightforward modifications, as long as this same expression is used in the actual Monte Carlo calculations for a double step in the imaginarytime direction. .
In order to match this spectrum, we take the derivative of Eq. (17) with respect to , evaluated between eigenstates of the proposed transfer matrix . We thus obtain, using the FeynmanHellmann theorem,
(30) 
where
(31) 
and
(32) 
The rest of the recipe consists in taking the righthand side of Eq. (30) and fitting the coefficients such that the first eigenvalues of the exact expression Eq. (29) (which correspond to the lowest eigenvalues prescribed by Lüscher’s formula) are reproduced. We may assume at this point that the coefficients are known, such that the eigenvectors are fixed when we set out to find the . In that case, the are determined by a linear system of equations of order :
(33) 
where
(34) 
and
(35) 
Once the coefficients have been determined, one can use Eq. (28) in a lattice calculation simply by replacing , and taking . The results of the fits for , are shown in Tables 1 and 2.
1  14.76869  –  –  –  – 

2  11.54894  1.74519  –  –  – 
3  10.74506  0.96946  0.40164  –  – 
4  10.31974  0.82605  0.17494  0.10404  – 
5  10.03874  0.74266  0.16064  0.02948  0.02878 
As a first illustration of the level of improvement that can be achieved for the energy, Fig. 2 shows the difference between the approximate spectrum with various levels of improvement and the exact spectrum , as a function of , through the quantity , where . As expected, with each new parameter a new eigenvalue is reproduced, with the concomitant reduction in the error.
Surprisingly, the approximations perform well considerably beyond the lowest eigenvalues that are fit. Indeed, with each new level of improvement we see, apart from a dramatic reduction in the error for the target eigenvalues, an extra reduction in the error that is evident even beyond .
iv.2 Contact
One of the many ways to define the contact (see Refs. ShinaContact ; BraatenReview ) is through the derivative of the energy with respect to the inverse scattering length:
(36) 
Since may be obtained directly from the logarithm of the partition function, we are again in a situation where we require a derivative of the transfer matrix with respect to a parameter, in this case . The generalization of the above definition to finite temperature in the grand canonical ensemble is
(37) 
where is the grand thermodynamic potential, is the temperature and is the chemical potential.
In manybody lattice calculations, using these definitions involves the following expression:
(38) 
where
(39) 
The first step towards using these expressions in combination with the improvement procedure is to take a formal derivative of with respect to in the twoparticle space, which we can treat exactly:
(40) 
In order to compute the change in the exact twoparticle energy due to a small change in the inverse scattering length, we use the fact that the energies are implicitly defined as solutions of Eq. (2), which implies
(41) 
where , and the derivative on the righthand side is to be evaluated at the corresponding solution of Eq. (2). Table 4 in the Appendix shows the first 30 roots of and the corresponding values of . In the derivation of Eq. (41), we have assumed that all the effectiverange parameters other than the scattering length are kept constant.
Having the exact target spectrum, we proceed by finding the corresponding expression in terms of the HS function , which we obtain using Eq. (17) and the FeynmanHellmann theorem:
(42) 
where, as before, we expand in terms of our chosen set of operators,
(43) 
and we determine the coefficients by fitting the diagonal matrix elements in the righthand side of Eq. (42) to the exact spectrum of Eqs. (40) and (41). As with the energy, the fitting procedure can be reduced to solving a set of linear equations of order . Illustrative results of such a fit are shown in Table 3.
1  0.36773  –  –  –  – 

2  0.14532  0.07568  –  –  – 
3  0.11370  0.02220  0.01957  –  – 
4  0.09659  0.01695  0.00415  0.00538  – 
5  0.08205  0.01278  0.00406  0.00023  0.00180 
As with the energy, a first glimpse at the level of improvement that can be obtained at this point. This is shown in Fig. 3, where we display the difference between the spectrum with various levels of improvement and the exact spectrum, divided by the latter. As in Fig. 2, each new parameter allows one to fit a new eigenvalue to high accuracy, matching the desired physics beyond the lowest momentum.
Unlike in Fig. 2, the improvement for eigenvalues beyond those explicitly fit is limited, breaking down after the eighth or ninth eigenvalue. From that point on toward higher energies, little improvement, if any, is observed as the order of the expansion is increased. This behavior is only unexpected in the light of Fig. 2, where the situation is (surprisingly) much more favorable. Possibly, part of the reason for this behavior is that the target spectra for the energy and the transfer matrix are monotonic (either increasing or decreasing) at low energies, whereas the one for the contact is closer to a random sequence (see Table 4 in the Appendix).
On the other hand, as mentioned before, there is no reason to believe that a given choice of operators will be equally useful for all observables. It is therefore natural to suspect that a different set of operators could yield a better expansion for the contact. This possibility remains to be studied.
iv.3 Extrapolation to the ground state
Taking a Slaterdeterminant state as a starting point, it is easy to see that the probability sum in groundstate calculations, for a temporal extent can be written as
(44) 
where are the exact energy eigenvalues and
(45) 
Taking a derivative of with respect to it is easy to see that in the large limit one can write
(46) 
where is the groundstate energy, , and
(47) 
Similarly, taking a derivative with respect to one can see that the first few leading contributions to the asymptotic behavior at large are given by
(48) 
where is the groundstate contact,
(49)  
(50) 
We shall use these expressions to motivate the extrapolations to the limit below.
V Illustrative results and conclusions
To illustrate the effect of improved operators in realistic lattice calculations, this section presents the results of groundstate MC calculations of the energy (Fig. 4) and the contact (Fig. 5), for an unpolarized system at unitarity.
Results are shown for 80 particles (40 per spin) in a volume of lattice points, for various levels of improvement . In each case, calculations were performed for time directions of extent (corresponding to roughly between 40 and 200), subsequently extrapolating to the limit. We have taken in lattice units, and a Slater determinant of plane waves as the starting guess for the groundstate wavefunction. For each value of , we obtained approximately samples of the auxiliary field . The statistics is enhanced by a factor of by the fact that the operators we have defined utilize every other time slice. Some obvious fluctuations remain, as evident from some degree of jaggedness in the data. This can be resolved by increasing the statistics, but it does not affect our main conclusions.
Once the improvements are turned on, i.e. for , we see that the change from to results in a considerable reduction in the energy. Our extrapolations to the ground state yield for the unimproved case, and for , respectively. This change, of roughly , is in line with our expectations based on the largescale calculations of Ref. Carlsonetal . While these results are still roughly greater than those of Ref. Carlsonetal , which found , it is remarkable that this can be achieved with a lattice. The remaining effects must be a combination of finite volume and finite imaginarytime step effects.
A similar situation is observed for the contact. Using the extrapolation formula of the previous section at leading plus nextto leading order, we obtain in the unimproved case, and and for , respectively. (We have discarded three data points at the lowest values of for the purpose of capturing the asymptotic behavior at large ; the fits are stable against further removal of points at low ). While the latter extrapolations clearly overlap when taking the uncertainty into account, they do not overlap with the unimproved case, which is clearly consistent with what we see in Fig. 5 at large . The remaining systematic errors, likely due to volume effects for the most part, appear to be only as large as , if we consider the most recent estimates of the contact in the ground state (, see Ref. Stefano ). As in the case of the energy, it is remarkable that such a small volume as already yields a result that is quite close to the best current estimate.
In conclusion, we have presented a methodology based on the work of Refs. Endresetal ; Endresetal2 , whereby one can generate not only improved transfer matrices but also improved operators, which account for finiterange effects in a systematic fashion. The improvement program can be carried out in a Galilean invariant or in a Galilean noninvariant way. We have chosen the latter because it provides a better approach to the unitary limit. As shown in Appendix B, the difference in the eigenvalues of the invariant and noninvariant transfer matrices is very small, such that the effect can be safely considered to be negligible. This is likely due to the fact that the lattice already breaks Galilean invariance. With the chosen improvements in place, we see a large change for the energy but a rather small change for the contact. Indeed, once finiterange effects are under control, the dominant contribution to systematic uncertainties is expected to be given by finitevolume effects. If so, the latter appear to be relatively small both for the energy (roughly ) as well as for the contact (roughly ).
While we have focused on the unitary limit, the tuning procedure for the transfer matrix and the operators can easily be applied to systems away from unitarity as well as away from the zero effectiverange limit. In combination with hybrid Monte Carlo (HMC) techniques, which have recently made it possible to access larger volumes than ever before (as large as ), we expect the method presented here to provide a powerful strategy to tackle the nonrelativistic manyfermion problem. This applies in particular to finitetemperature calculations where lowering the ultraviolet cutoff (while keeping effectiverange effects under control) reduces the size of the basis, speeding up the computations considerably.
In this regard, it should also be pointed out that improving the transfer matrix affects the performance of the HMC algorithm only marginally. Indeed, the scaling with system size (particle number, spacetime volume) remains unchanged when increasing . Only the prefactor in the scaling law increases somewhat (by about ) due to the need for extra Fourier transforms when applying the operator . Aside from this, the force in the molecular dynamics (MD) part of the HMC algorithm becomes somewhat larger when is improved, such that somewhat smaller MD time steps (not to be confused with ) are needed. In short, the performance of the HMC algorithm is largely unaffected by the use of improved transfer matrices.
Finally, we would like to stress that the operator proposed here for the calculation of the contact (based on the derivative of the twobody energy, in turn computed using Eq. (41)) constitutes a novel way to calculate in MC calculations. This method should be contrasted with others based on extracting the derivative of the equation of state via finite differences, using the momentum distribution, or using the pair distribution function. Each of these will, in general, behave differently in terms of their systematic effects.
Acknowledgements.
Thanks are in order to J. Carlson, W. Detmold, T. Lähde, D. Lee, Y. Nishida, E. Passemar, and S. Tan for discussions and/or critical comments on the manuscript. I am especially indebted to E. R. Anderson, for numerous discussions that led to this work, and to A. N. Nicholson for clarifications on the work of Ref. Endresetal . I would also like to thank the Institute for Nuclear Theory in Seattle, WA, as well as the “Quarks, Hadrons and Nuclei” group at the University of Maryland, where some of the early work and discussions took place. This work was supported in part by U.S. Department of Energy, Office of Nuclear Physics under contract DEFC0207ER41457 (UNEDF SciDAC).Appendix A Evaluation of .
There are many ways to evaluate the function efficiently. The one that I have found easiest to implement and understand is essentially the one communicated to me by Shina Tan, which I reproduce here, with some minor modifications that I have found useful.
The first step is to enhance the convergence of the sum by introducing an exponential factor and separating the singularity in the tail (see Ref. Shina ):
(51) 
where , the sums are over all triplets of integers , and is a small positive parameter. The second sum has no convergence problems, in fact it converges very quickly, so we shall put it aside by defining
(52) 
1  0.0959007  123.82387 

2  0.4728943  39.75514 
3  1.4415913  82.36519 
4  2.6270076  106.24712 
5  3.5366199  84.23133 
6  4.2517060  161.88763 
7  5.5377008  212.49220 
8  7.1962632  62.95336 
9  8.2879537  231.79580 
10  9.5345314  247.82611 
11  10.5505341  233.82976 
12  11.7014957  185.61411 
13  12.3102392  183.65019 
14  13.3831152  316.68684 
15  15.3537375  82.86757 
16  16.1218253  506.59914 
17  17.5325415  371.40245 
18  18.6053932  308.00372 
19  19.5186394  255.97969 
20  20.4033187  329.98905 
21  21.6944179  394.81924 
22  23.0194727  94.98929 
23  24.3306210  342.25749 
24  25.3016129  526.27127 
25  26.6803600  514.90705 
26  27.8780019  150.20773 
27  29.6156511  548.38017 
28  31.3536974  114.02114 
29  32.1958982  443.21169 
30  33.4483351  452.78989 

The first sum, on the other hand, converges just as slowly as the original one (once we subtract the term). Let us then analyze this first sum. For , we may approximate it very accurately in terms of an integral:
(53) 
where we have used the spherical symmetry of the integrand. Writing
(54) 
we separate the integrals into two terms. First,
(55) 
where the cancels with the term in Eq. (3), and the integral in the second term is equal to in the limit . Second, we have
(56) 
To treat this term we write
(57) 
and exchange the order of integration of and , to obtain
(58) 
all of which can be easily evaluated in the limit , where the full result can be written as
(59) 
which is to be evaluated for .
With the above expressions at hand, it is straightforward to find the first derivative of with respect to :
(60) 
where
(61) 
Table 4 shows the first 30 roots of , along with the corresponding values of the derivative .
Appendix B Galilean invariance vs. noninvariance.
One cannot fail to notice that not only is the transfer matrix of Eq. (16) not symmetric, an issue we managed to deal with above, it is also not fully Galilean invariant. Indeed, while translation and rotation invariance are preserved (in their discrete forms set by the lattice), Galileanboost symmetry is broken. This is because we have assumed that was promoted to an operator that acts on everything that appears to its right.
On the other hand, if we assume that acts on the auxiliary field only, as in Ref. Endresetal , one arrives directly at a local, Galileaninvariant lattice theory. Indeed, in that case Eq. (16) changes form in a very simple manner, namely the function is replaced by , and similarly for the other spin. The function has then a clear physical significance as the momentum transfer in a twobody collision.
Apart from the obvious desirability of Galilean invariance, these observations seem to indicate that one should choose this form over the noninvariant version. However, when comparing the approach to the unitary point based on improvements, as shown in Figs. 1 and 6, it is clear that the noninvariant version is much superior to its invariant counterpart. Furthermore, comparing with the results of Ref. Endresetal , it seems clear that one would need much larger lattices in the Galilean invariant version to make Fig. 6 look like Fig. 1, which we have checked with our codes. (Note that Fig. 2 in Ref. Endresetal presents the same kind of plot but corresponding to a much larger lattice size, namely .)
The reason for this difference is likely due to the fact that (discrete) Galilean invariance is only respected in an infinite lattice volume. Indeed, in any finite lattice the centerofmass and relative motions are actually not separable. Boosting to frames of different total momentum will, in general, result in Hilbert spaces of vastly different dimensions for the subspace of relative motion. In particular, the zero total momentum Hilbert space, used to tune our twobody interaction, has by far the highest dimension. Based on this assessment, we decided to focus on the noninvariant form in this work.
Should one choose to implement the Galilean invariant form instead, as in Ref. Endresetal , there are two technical details that should not be overlooked. First, in Eq. (16) the function becomes and therefore needs to be evaluated outside the singleparticle momentum lattice where we originally defined it (c.f. Eq. (13)). The extension to the larger domain should be done respecting the periodicity of the momentum lattice.
Second, if such a periodicity is to be reconciled with the continuum form of the noninteracting dispersion relation, which we assumed to be , we should impose a spherically symmetric cutoff in momentum space. This last condition should have little impact on the results in practice, as long as the systems under consideration are somewhat dilute. From the point of view of computational performance, however, such a restriction on phase space results in important gains, particularly at finite temperature, where all the elements in the basis are evolved in imaginary time.
Regardless of the form of the action, it is possible to evaluate the degree to which Galilean invariance is broken by taking the improved transfer matrix and computing its eigenvalues in a moving frame. As long as Galilean invariance is respected, inserting the eigenvalues into the RummukainenGottlieb formula RummuGott should yield the same scattering phase shift as Lüscher’s formula. This will of course include breaking effects coming from both the action and the lattice itself.
Here, we limit ourselves to comparing the eigenvalues obtained from the invariant and noninvariant formulations, in various frames. In Fig. 7 we show a plot of the relative eigenvalue difference between the Galilean invariant and noninvariant improved transfer matrices, evaluated in various moving frames (i.e. nonzero centerofmass momentum). As can be appreciated in that plot, the difference between the invariant and noninvariant improvements is extremely small, which shows that one can use the noninvariant version of the improvement program without introducing a noticeable systematic effect.
References
 (1) I. Montvay, G. Münster, Quantum fields on the lattice (Cambridge University Press, New York, NY, 1994); H. Rothe, Lattice Gauge Theories, 3rd ed. (World Scienti c, Singapore, 2005). T. DeGrand and C. DeTar, Lattice Methods for Quantum Chromodynamics (World Scienti c, Singapore, 2006).
 (2) P. Lepage, Lectures given at the VIII Jorge Andre Swieca Summer School (Brazil, 1997); arXiv:nuclth/9706029. G. P. Lepage, in Proceedings of TASI’89: From Actions to Answers, T. DeGrand and D. Toussaint (eds.) (World Scientific, 1989); arXiv:hepph/0506330v1.
 (3) K. Symanzik, Nucl. Phys. B 226, 187 (1983); ibid. 226, 205 (1983).
 (4) M. G. Endres, D. B. Kaplan, J.W. Lee, A. N. Nicholson, Phys. Rev. A 84, 043644 (2011).
 (5) M. G. Endres, D. B. Kaplan, J.W. Lee, A. N. Nicholson, PoS (Lattice 2010) 182. J.W. Lee, M. G. Endres, D. B. Kaplan, A. N. Nicholson, PoS (Lattice2010) 197. A. N. Nicholson, M. G. Endres, D. B. Kaplan, J.W. Lee, PoS (Lattice2010) 206.
 (6) Ultracold Fermi Gases, Proceedings of the International School of Physics “Enrico Fermi”, Course CLXIV, Varenna, June 20 – 30, 2006, M. Inguscio, W. Ketterle, C. Salomon (Eds.) (IOS Press, Amsterdam, 2008).
 (7) M. Lüscher, Commun. Math. Phys. 105, 153 (1986).
 (8) R. L. Stratonovich, Sov. Phys. Dokl. 2 (1958) 416; J. Hubbard, Phys. Rev. Lett. 3 (1959) 77.
 (9) J.W. Negele and H. Orland, Quantum ManyParticle Systems (AddisonWesley, Redwood City, CA, 1988).
 (10) A. Privitera, M. Capone, C. Castellani, Phys. Rev. B 81, 014523 (2010); A. Privitera, M. Capone, Phys. Rev. A 85, 013640 (2012).
 (11) J. Carlson, S. Gandolfi, K. E. Schmidt, S. Zhang, Phys. Rev. A 84, 061602(R) (2011).
 (12) S. Bour, X. Li, D. Lee, UlfG. Meissner, and L. Mitas, Phys. Rev. A 83, 063619 (2011).
 (13) J. Braun, S. Diehl, and M. M. Scherer, Phys. Rev. A 84, 063616 (2011).
 (14) A. Bazabov, D. Toussaint, C. Bernard, J. Laiho, C. DeTar, L. Levkova, M. B. Oktay, U. M. Heller, J. E. Hetrick, P. B. Mackenzie, R. Sugar, R. S. Van de Water, Rev. Mod. Phys. 82, 1349 (2010).
 (15) A. Bulgac, J. E. Drut, P. Magierski, Phys. Rev. Lett. 99 (2007) 120401; Phys. Rev. Lett. 96 (2006) 090404; Phys. Rev. A 78 (2008) 023625; P. Magierski, G. Wlazłowski, A. Bulgac, J. E. Drut, Phys. Rev. Lett. 103 (2009) 210403.
 (16) J. E. Drut, T. A. Lähde, T. Ten, Phys. Rev. Lett. 106, 205302 (2011).
 (17) J. E. Drut, T. A. Lähde, G. Wlazłowski,P. Magierski, Phys. Rev. A 85, 051601(R) (2012).
 (18) R. Blankenbecler, D. J. Scalapino, R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
 (19) D. Lee, Phys. Rev. C 78, 024001 (2008); Prog. Part. Nucl. Phys. 63, 117 (2009).
 (20) S. Tan, Ann. Phys. 323, 2952 (2008); ibid. 323, 2971 (2008); ibid. 323, 2987 (2008).
 (21) E. Braaten, Universal Relations for Fermions with Large Scattering Length, in The BCSBEC crossover and the Unitary Fermi Gas W. Zwerger (ed.) (Springer, Berlin, 2011); arXiv:1008.2922.
 (22) S. Tan, Phys. Rev. A 78, 013636 (2008).
 (23) J. ZinnJustin, Quantum Field Theory and Critical Phenomena (Oxford University Press, New York, 2002).
 (24) S. Gandolfi, K. E. Schmidt, J. Carlson, Phys. Rev. A 83, 041601(R) (2011).
 (25) K. Rummukainen, S. Gottlieb, Nucl. Phys. B 450, 397 (1995); C.h. Kim, C.T. Sachrajda, S. R. Sharpe, Nucl. Phys. B 727, 218 (2005).