Lattice methods and effective field theory

Lattice methods and effective field theory

Amy N. Nicholson Department of Physics, University of California, Berkeley, Berkeley CA 94720, USA
August 22, 2019

Lattice field theory is a non-perturbative tool for studying properties of strongly interacting field theories, which is particularly amenable to numerical calculations and has quantifiable systematic errors. In these lectures we apply these techniques to nuclear Effective Field Theory (EFT), a non-relativistic theory for nuclei involving the nucleons as the basic degrees of freedom. The lattice formulation of Endres et al. (2011a, 2013) for so-called pionless EFT is discussed in detail, with portions of code included to aid the reader in code development. Systematic and statistical uncertainties of these methods are discussed at length, and extensions beyond pionless EFT are introduced in the final Section.

I Introduction

Quantitative understanding of nuclear physics at low energies from first principles remains one of the most challenging programs in contemporary theoretical physics research. While physicists have for decades used models combined with powerful numerical techniques to successfully reproduce known nuclear structure data and make new predictions, currently the only tools available for tackling this problem that have direct connections to the underlying theory, Quantum Chromodynamics (QCD), as well as quantifiable systematic errors, are Lattice QCD and Effective Field Theory (EFT). In principle, when combined these techniques may be used to not only quantify any bias introduced when altering QCD in order to make it computationally tractable, but also to better understand the connection between QCD and nuclear physics.

The lattice is a tool for discretizing a field theory in order to reduce the path integral, having an infinite number of degrees of freedom, to a finite-dimensional ordinary integral. After rendering the dimension finite (though extremely large), the integral may then be estimated on a computer using Monte Carlo methods. Errors introduced through discretization and truncation of the region of spacetime sampled are controlled through the spatial and temporal lattice spacings, , and the number of spatial and temporal points, . Thus, these errors may be quantified through the lattice spacing dependence of the observables, and often may be removed through extrapolation to the continuum and infinite volume limits.

LQCD is a powerful and advanced tool for directly calculating low-energy properties of QCD. However, severe computational issues exist when calculating properties of systems with nucleons. Unfortunately, these problems grow rapidly with the number of nucleons in the system.

The first issue is the large number of degrees of freedom involved when using quark fields to create nucleons. In order to calculate a correlation function for a single nucleon in LQCD using quarks (each of which has twelve internal degrees of freedom given by spin and color), one has to perform all possible Wick contractions of the fields in order to build in fermion antisymmetrization. For example, to create a proton using three valence quark operators requires the calculation of two different terms corresponding to interchanging the two up quark sources. The number of contractions involved for a nuclear correlation function grows with atomic number and mass number as . For He this corresponds to terms111This is a very naïve estimate; far more sophisticated algorithms exist with power-law scaling.!

The second major problem occurs when performing a stochastic estimate of the path integral. A single quark propagator calculated on a given gauge field configuration may be a part of either a light meson or a heavy nucleon. However, the difference cannot be determined until correlations with the other quark fields present are built in by summing over a sufficiently large number of these field configurations222This interpretation of the signal-to-noise problem has been provided by David B. Kaplan.. This leads to large fluctuations from configuration to configuration, and a stochastic signal-to-noise ratio, , which degrades exponentially with the number of nucleons in the system,


where is the nucleon mass and is the pion mass Lepage (1989). This is currently the major limiting factor for the size of nuclear which can be probed using LQCD. The best calculations we have from LQCD using multiple nucleons to date are in the two-nucleon sector Berkowitz et al. (2015); Kurth et al. (2015); Nicholson et al. (2015); Orginos et al. (2015); Detmold et al. (2016); Chang et al. (2015); Beane et al. (2015a, b, 2014, 2013a, 2013b, 2012a, 2010); Yamazaki (2015); Yamazaki et al. (2015, 2014); Yamazaki et al. (2012a, b); Doi et al. (2015a, b); Ishii et al. (2007); Murano et al. (2014); Aoki (2013); Murano (2013); Ishii et al. (2012); Inoue et al. (2010), while fewer calculations have been performed for three and four nucleon systems Beane et al. (2009a, 2013b, 2014, 2015b); Chang et al. (2015); Yamazaki (2015); Yamazaki et al. (2015, 2014); Yamazaki et al. (2012a, b); Doi et al. (2012); however, even for two nucleon systems unphysically large pion masses must be used in order to reduce the noise problem. We will discuss signal-to-noise problems in more detail in Sec. III.1.

Starting from an EFT using nucleons as the fundamental degrees of freedom greatly reduces the consequences from both of these issues. EFTs also enjoy the same benefit as the lattice over traditional model techniques of having quantifiable systematic errors, this time controlled by the cutoff of the EFT compared to the energy regime studied. For chiral EFTs this scale is generally MeV. Systematic errors can be reduced by going to higher orders in an expansion of , where is the momentum scale probed, with the remaining error given by the size of the first order which is not included. In a potential model there is no controlled expansion, and it is generally unknown how much the results will be affected by leaving out any given operator. In addition, field theories provide a rigorous mathematical framework for calculating physical processes, and can be directly translated into a lattice scheme.

In these lecture notes we will explore the use of lattice methods for calculating properties of many-body systems starting from nuclear EFT, rather than QCD. Our discussion will begin with understanding a very basic nuclear EFT, pionless EFT, at leading order. We will then proceed to discretize this theory and set up a framework for performing Monte Carlo calculations of our lattice theory. We will then discuss how to calculate observables using the lattice theory, and how to understand their associated statistical uncertainties. Next we will discuss quantifying and reducing systematic errors. Then we will begin to add terms to our theory going beyond leading order pionless EFT. Finally, we will discuss remaining issues and highlight some successes of the application of these methods by several different groups.

Ii Basics of Effective Field Theory and Lattice Effective Field Theory

ii.1 Pionless Effective Field Theory

To develop an EFT we will first write down all possible operators involving the relevant degrees of freedom within some energy range (determined by the cutoff) that are consistent with the symmetries of the underlying theory. Each operator will be multiplied by an unknown low-energy constant which may be fixed by comparing an observable with experiment or lattice QCD. In order to reduce this, in principle, infinite number of operators to a finite number we must also establish a power-counting rule for neglecting operators that do not contribute within some desired accuracy. This is a notoriously difficult problem for nuclear physics, and is in general observable and renormalization scheme dependent. Here, we will only briefly touch upon two common power-counting schemes, the so-called Weinberg and KSW expansions Weinberg (1990, 1991); Kaplan et al. (1996, 1998a, 1998b). For reviews of these and other power-counting schemes, see Epelbaum et al. (2009); Epelbaum (2010); Machleidt and Entem (2011).

The simplest possible nuclear EFT involves non-relativistic nucleon fields interacting via delta functions. This is known as a pionless EFT, and is only relevant for energy scales up to a cutoff . Below this scale, the finite range of pion exchange cannot be resolved, and all interactions appear to be point-like. In this discussion we will closely follow that of Ref. Kaplan (2005). For the moment, let’s just consider a theory of two-component (spin up/down) fermion fields, , with the following Lagrangian,




represents the nucleon mass, are unknown, low-energy constants (LECs) which may be fixed by comparing to experimental or LQCD results, and all spin indices are suppressed. Because the effective theory involves dynamical degrees of freedom that are only relevant up to a certain scale, we must define a cutoff, , above which the theory breaks down. In general, the LECs scale as , where dim represents the dimension of the operator associated with the LEC. According to naïve power counting, the term in Eq. (2) should be suppressed relative to the term, because adding a derivative to an operator increases its dimension. One should be careful in practice, however, because naïve power counting does not always hold, as we will see several times throughout these lectures.

ii.1.1 Two particle scattering amplitude

In order to set the coefficients , we may look to experimental scattering data. In particular, if we wish to set the coefficient we should consider two-particle -wave scattering because the operator associated with contains no derivatives. and other LECs may be set using - and higher-wave scattering data. Recall that the S-matrix for non-relativistic scattering takes the following form:


where is the scattering momentum and is the scattering amplitude. For -wave scattering the amplitude may be written as,


where is the -wave scattering phase shift. Given a short-range two-body potential, the scattering phase shift has a well-known expansion for low momenta, called the effective range expansion,


where is the scattering length, is the effective range, and and higher order terms are referred to as shape parameters. The effective range and shape parameters describe the short-range details of the potential, and are generally of order of the appropriate power of the cutoff in a naturally tuned scenario.

The scattering length may be used to describe the asymptotic behavior of the radial wavefunction. In particular, consider two-particles interacting via an attractive square-well potential. If the square-well is sufficiently strongly attractive, the wavefunction turns over and goes to zero at some finite characteristic length. This means the system is bound and the size of the bound state is given by the scattering length, . On the other hand, if the wavefunction extends over infinite space, then the system is in a scattering state and the scattering length may be determined as the distance from the origin where the asymptote of the wavefunction intersects the horizontal axis (see Fig. 1). This implies that the scattering length in the case of a scattering state is negative. If the potential is tuned to give a system which is arbitrarily close to the crossover point from a bound state to a scattering state, corresponding to infinite scattering length, the state is described as being near unitarity, because the unitarity bound on the scattering cross section is saturated at this point. Note that this implies that the scattering length may be any size and is not necessarily associated with the scale set by the cutoff. However, such a scenario requires fine-tuning of the potential. Such fine-tuning is well-known to occur in nuclear physics, with the deuteron and neutron-neutron -wave scattering being notable examples.

Figure 1: Sketches of two-body radial wavefunctions vs. corresponding to various scattering lengths. From left to right: , ,.

A many-body system composed of two-component fermions with an attractive interaction is known to undergo pairing between the species (higher -body interactions are prohibited by the Pauli exclusion principle), such as in neutron matter, found in the cores of neutron stars, which is composed of spin up and spin down neutrons. At low temperature, these bosonic pairs condense into a coherent state. If the interaction is only weakly attractive, the system will form a BCS state composed of widely separated Cooper pairs, where the average pair size is much larger than the average interparticle spacing. On the other hand, if the interaction is strongly attractive then the pairs form bosonic bound states which condense into a Bose-Einstein condensate. The crossover between these two states corresponds to the unitary regime, and has been studied extensively in ultracold atom experiments, where the interaction between atoms may be tuned using a Feshbach resonance. In this regime, the average pair size is equal to the interparticle spacing (given by the inverse density), which defines the only scale for the system. Thus, all dimensionful observables one wishes to calculate for this system are determined by the appropriate power of the density times some dimensionless constant. For a review of fermions in the unitary regime, see e.g., Giorgini et al. (2008); Block et al. (2008).

ii.1.2 Two-body LECs

Returning to our task of setting the couplings using scattering parameters as input, we might consider comparing Eq. (2) and Eq. (6), to determine the LEC using the scattering length, using the effective range, and so forth. To see how this is done in practice we may compute the scattering amplitude in the effective theory, and match the coefficients to the effective range expansion. Let’s begin using only the first interaction term in the effective theory, corresponding to . Diagrammatically, the scattering amplitude may be written as the sum of all possible bubble diagrams (see Fig. 2). Because the scattering length may take on any value, as mentioned previously, we cannot assume that the coupling is small, so we should sum all diagrams non-perturbatively. The first diagram in the sum is given by the tree level result, . If we assume that the system carries energy , then the second diagram may be labeled as in Fig. 3, and gives rise to the loop integral,


Performing the integral over and the solid angle gives


where I have introduced a hard momentum cutoff, . Removing the cutoff by taking it to infinity results in


Because the interaction is separable, the th bubble diagram is given by products of this loop function. Thus, the scattering amplitude is factorizable, and may be written


We may now compare Eqs. (5,6) and Eq. (11) to relate the coupling to the scattering phase shift. This is easiest to do by equating the inverse scattering amplitudes,


where I have used Eq. (6) cut off at leading order. We now have the relation


between the coupling and the physical scattering length.

Figure 2: Two-body scattering amplitude represented as a sum of bubble diagrams corresponding to a single contact interaction with coupling .
Figure 3: Feynman diagram for a single bubble in Fig. 2, giving rise to the loop integral Eq. (7).

Note that the coupling runs with the scale ; the particular dependence is determined by the regularization and renormalization scheme chosen. In order to understand the running of the coupling we may examine the beta function. To do so we first define a dimensionless coupling,


then calculate


This function is a simple quadratic that is plotted in Fig. 4. The beta function has two zeroes, , corresponding to fixed points of the theory. At a fixed point, the coupling no longer runs with the scale , and the theory is said to be scale-invariant (or conformal, given some additional conditions). This means that there is no intrinsic scale associated with the theory. The fixed point at is a trivial fixed point, and corresponds to a non-interacting, free field theory (zero scattering length). The other, non-trivial fixed point at corresponds to a strongly interacting theory with infinite scattering length; this is the unitary regime mentioned previously. Here, not only does the scattering length go to infinity, as does the size of the radial wavefunction, but the energy of the bound state (as approached from ) goes to zero and all relevant scales have vanished. Note that this is an unstable fixed point; the potential must be finely tuned to this point or else the theory flows away from unitarity as (IR limit).

Figure 4: Beta function (Eq. (16)) for the two-body contact interaction. Arrows represent the direction of flow toward the IR.

Generally perturbation theory is an expansion around free field theory, corresponding to a weak coupling expansion. This is the approach used as part of the Weinberg power counting scheme for nuclear EFT Weinberg (1990, 1991). However, in some scattering channels of interest for nuclear theory the scattering length is indeed anomalously large, such as the and nucleon-nucleon scattering channels, where


Such large scattering lengths suggest that an expansion around the strongly coupled fixed-point of unitarity may be a better starting point and lead to better convergence. This approach was taken by Kaplan, Savage, and Wise and led to the KSW power-counting scheme Kaplan et al. (1998b, a, 1996). Unfortunately, nuclear physics consists of many scales of different sizes and a consistent power-counting framework with good convergence for all observables has yet to be developed; in general the convergence of a given scheme depends on the scattering channels involved.

Because nuclear physics is not weakly coupled in all channels, non-perturbative methods, such as lattice formulations, will be favorable for studying few- and many-body systems, where two-body pairs may interact through any combination of channels simultaneously. Due to the scale-invariant nature of the unitary regime, it provides a far simpler testbed for numerical calculations of strongly-interacting theories, so we will often use it as our starting point for understanding lattice EFT methods.

ii.2 Lattice Effective Field Theory

Our starting point for building a lattice EFT will be the path integral formulation of quantum field theory in Euclidean spacetime. The use of Euclidean time allows the exponent of the path integral to be real (in certain cases), a property which will be essential to our later use of stochastic methods for its evaluation. Given a general theory for particles obeying a Lagrangian density


where is the Euclidean time, the chemical potential, and is the Hamiltonian density, the Euclidean path integral is given by


If the integral over Euclidean time is compact, then the finite time extent acts as an inverse temperature, and we may draw an analogy with the partition function in statistical mechanics, . This analogy is often useful when discussing lattice formulations of the path integral. In this work we will generally consider and create non-zero particle density by introducing sources and sinks for particles and calculating correlation functions.

We discretize this theory on a square lattice consisting of points, where is the number of points in all spatial directions, and is the number of temporal points. We will focus on zero temperature physics, corresponding to large 333The explicit condition on required for extracting zero temperature observables will be discussed in Sec. III. We must also define the physical distance between points, the lattice spacings , where by dimensional analysis for non-relativistic theories. The fields are now labeled by discrete points, , and continuous integrals are replaced by discrete sums, .

ii.2.1 Free field theory

To discretize a free field theory, we must discuss discretization of derivatives. The simplest operator which behaves as a single derivative in the continuum limit is a finite difference operator,


where is a unit vector in the -direction. The discretized second derivative operator must involve two hops, and should be a symmetric operator to behave like the Laplacian. A simple possibility is


We can check the continuum limit by inspecting the corresponding kinetic term in the action,


The fields may be expanded in a plane wave basis,


for spatial indices, , leading to


After performing the sum over the first piece in brackets gives , while the second is proportional to , resulting in,


Finally, expanding the sine function for small gives,


where I’ve used the finite volume momentum to rewrite the expression in square brackets. Thus, we have the correct continuum limit for the kinetic operator. Note that for larger momenta, approaching the continuum limit requires smaller . However, this is only one possibility for a kinetic term. We can always add higher dimension operators (terms with powers of in front of them), in order to cancel leading order terms in the expansion Eq. (28). This is a form of what’s called improvement of the action, and will be discussed in more detail in Sec. IV.

Adding a temporal derivative term,


we can now write down a simple action for a non-relativistic free-field theory,


where I’ve defined a matrix whose entries are blocks,


where contains the spatial Laplacian, and therefore connects fields on the same time slice (corresponding to diagonal entries of the matrix ), while the temporal derivative contributes the off-diagonal pieces. Note that the choice of “1” in the lower left corner corresponds to anti-periodic boundary conditions, appropriate for fermionic fields. For zero temperature calculations the temporal boundary conditions are irrelevant, and it will often be useful to choose different temporal boundary conditions for computational or theoretical ease.

ii.2.2 Interactions

Now let’s discuss adding interactions to the theory. We’ll focus on the first term in a nuclear EFT expansion, the four-fermion interaction:


where now explicitly label the particles’ spins (or alternatively, flavors). Because anti-commuting fields cannot easily be accommodated on a computer, they must be integrated out analytically. The only Grassmann integral we know how to perform analytically is a Gaussian, so the action must be bilinear in the fields. One trick for doing this is called a Hubbard-Stratonovich (HS) transformation, in which auxiliary fields are introduced to mediate the interaction. The key is to use the identity,


where I have dropped the spacetime indices for brevity. This identity may be verified by completing the square in the exponent on the right hand side and performing the Gaussian integral over the auxiliary field . This form of HS transformation has the auxiliary field acting in what is called the density channel . It is also possible to choose the so-called BCS channel, , the usual formulation used in BCS models, however this causes a so-called sign problem when performing Monte Carlo sampling, as will be discussed in detail in Sec. III.1.1. Transformations involving non-Gaussian auxiliary fields may also be used, such as

compact continuous: (35)

These formulations may have different pros and cons in terms of computational and theoretical ease for a given problem, and should be chosen accordingly. For example, the interaction is conceptually and computationally the simplest interaction, however, it also induces explicit and higher-body interactions in systems involving more than two-components which may not be desired.

ii.2.3 Importance sampling

The action may now be written with both kinetic and interaction terms,


where the matrix includes blocks which depend on the auxiliary field , and also contains non-trivial spin structure that has been suppressed. The partition function can be written


where the integration measure for the field, , depends on the formulation chosen,


With the action in the bilinear form of Eq. (36), the fields can be integrated out analytically, resulting in


Observables take the form


Through the use of discretization and a finite volume, the path integral has been converted into a standard integral with finite dimension. However, the dimension is still much too large to imagine calculating it on any conceivable computer, so we must resort to Monte Carlo methods for approximation. The basic idea is to generate a finite set of field configurations of size according to the probability measure , calculate the observable on each of these configurations, then take the mean as an approximation of the full integral,


Assuming the central limit theorem holds, for large enough (a non-trivial condition, as will be discussed in Sec. III.2), the distribution of the mean approaches a Gaussian, and the error on the mean falls off with the square root of the sample size.

There are several algorithms on the market for generating field configurations according to a given probability distribution, and I will only briefly mention a few. Lattice calculations are particularly tricky due to the presence of the determinant in Eq. (39), which is a highly non-local object and is very costly to compute. One possible algorithm to deal with this is called determinantal Monte Carlo, which implements local changes in , followed by a simple Metropolis accept/reject step. This process can be rather inefficient due to the local updates. An alternative possibility is Hybrid Monte Carlo, commonly used for lattice QCD calculations, in which global updates of the field are produced using molecular dynamics as a guiding principle. Note that the field must be continuous in order to use this algorithm due to the use of classical differential equations when generating changes in the field. Also common in lattice QCD calculations is the use of pseudofermion fields as a means for estimating the fermion determinant. Here the determinant is rewritten in terms of a Gaussian integral over bosonic fields, ,


This integral is then evaluated stochastically. These are just a sample of the available algorithms. For more details on these and others in the context of non-relativistic lattice field theory, see Drut and Nicholson (2013).

ii.2.4 Example formulation

Now that we have developed a general framework for lattice EFT, let’s be explicit and make a few choices in order to further our understanding and make calculations simpler. The first choice I’m going to make is to use a field, so that is trivial. The next simplification I’m going to make is to allow the fields to live only on temporal links,


Note that we are free to make this choice, so long as the proper four-fermion interaction is regained in the continuum limit. This choice renders the interaction separable, as it was in our continuum effective theory. This means we may analytically sum two-body bubble chain diagrams as we did previously in order to set the coupling using some physical observable (see Fig. 5).

With this choice we can now write the -matrix explicitly as


where . Now the -dependence exists only on the upper diagonal, as well as the lower left due to the boundary condition. This block will be eliminated through our final choice: open boundary conditions in time for the fields, . As mentioned previously, we are free to choose the temporal boundary conditions as we please, so long as we only consider zero temperature (and zero chemical potential) observables.

Figure 5: Two-body scattering amplitude of Fig. 2, where the contact interaction has been replaced in the second line by exchange of a dimer auxiliary field via a Hubbard-Stratonovich transformation.

With this set of choices the matrix consists purely of diagonal elements, , and upper diagonal elements, . One property of such a matrix is that the determinant, which is part of the probability distribution, is simply the product of diagonal elements, . Note that is completely independent of the field . This means that the determinant in this formulation has no impact on the probability distribution , and therefore never needs to be explicitly computed, greatly reducing the computational burden. Thus in all of our calculations, performing the path integral over simply amounts to summing over at each lattice site.

Finally, this form of also makes the calculation of propagators very simple. The propagator from time 0 to may be written,


where , and all entries are matrices which may be projected onto the desired state. This form suggests a simple iterative approach to calculating propagators: start with a source (a spatial vector projecting onto some desired quantum numbers and interpolating wavefunction), hit it with the kinetic energy operator corresponding to free propagation on the time slice, then hit it with the field operator on the next time link, then another free kinetic energy operator, and so on, finally projecting onto a chosen sink vector.

As will be discussed further in Secsystematic, it is often preferable to calculate the kinetic energy operator in momentum space, while the auxiliary field in must be generated in position space. Thus, Fast Fourier Transforms (FFTs) may be used between each operation to quickly translate between the bases. Example code for generating source vectors, kinetic operators, and interaction operators will be provided in later Sections.

A cartoon of this process on the lattice is shown in Fig. 6. The choice of auxiliary fields also simplifies the understanding of how four-fermion interactions are generated. On every time link, imagine performing the sum over . If there is only a single fermion propagator on a given link this gives zero contribution because the term is proportional to . However, on time slices where two propagators overlap, we have instead . In sum, anywhere two fermions exist at the same spacetime point a factor of contributes, corresponding to an interaction.

Figure 6: Schematic of a lattice calculation for a two-particle correlation function. The two particles (red and blue lines) propagate through the lattice between source and sink , seeing particular values of the auxiliary field, , on each time link. If two particles occupy the same temporal link, then upon summation over all possible values of at each link, a non-zero contribution is generated by the interaction term because .

ii.2.5 Tuning the two-body interaction

There are several ways to set the two-body coupling. Here we will explore two methods, using different two-body observables. The first involves calculating the two-particle scattering amplitude, and tuning the coupling to reproduce known scattering parameters, to make a connection with our previous calculation for the effective theory. The second method uses instead the energy spectrum of a two-particle system in a box. This powerful method will be useful later when we begin to improve the theory in order to reduce systematic errors.

We have calculated the scattering amplitude previously for our effective theory using a momentum cutoff. For the first method for tuning the coupling, we will calculate it again using our lattice theory with the lattice cutoff as a regulator. First we need the single particle free propagator:


where I’ve set (we will use this convention from now on until we begin to discuss systematic errors), and have used the previously defined discretized Laplacian operator. I’ve written the propagator in a mixed representation, as this is often useful in lattice calculations for calculating correlation functions in time when the kinetic operator, , is diagonal in momentum space.

The diagrammatic two-particle scattering amplitude is shown on the bottom line in Fig. 5. Because we have chosen the interaction to be separable, the amplitude can be factorized:


where the one loop integral, , will be defined below. As before, in order to set a single coupling we need one observable, so we use the effective range expansion for the scattering phase shift to leading order,


Relating Eqs. (50,51), we find


We will now evaluate the loop integral using the free single particle propagators, Eq. (47),


This final sum may be calculated numerically for a given and (governing the values of momenta included in the sum), as well as for different possible definitions of the derivative operators contained in , giving the desired coupling, , via Eq. (52).

The second method for setting the coupling utilizes the calculation of the ground state energy of two particles. We start with the two-particle correlation function,


where is a source (sink) wavefunction involving one spin up and one spin down particle. Integrating out the fermion fields gives,


I will now write out the components of the matrices explicitly:


The first (last) piece in angle brackets represents the position space wavefunction created by the sink (source). All fields in Eq. (62) are uncorrelated, so we can perform the sum for each time slice independently. One such sum is given by,


where the cross terms vanish upon performing the sum. If we make the following definitions,


then we can write the two-particle correlation function as,


where I have made the definition


Recall from statistical mechanics that correlation functions may be written as insertions of the transfer matrix, , acting between two states,


Then we may identify in Eq. (68) as the transfer matrix of the theory, . This in turn implies that the logarithm of the eigenvalues of give the energies of the two-particle system.

We will now evaluate the transfer matrix in momentum space:


where I have made the definition,


The eigenvalues of the matrix may be evaluated numerically to reproduce the entire two-particle spectrum. However, for the moment we only need to set a single coupling, , so one eigenvalue will be sufficient. The largest eigenvalue of the transfer matrix, corresponding to the ground state, may be found using a simple variational analysis444Many thanks to Michael Endres for the following variational argument.. Choosing a simple trial state wavefunction,


subject to the normalization constraint,


we now need to maximize the following functional:


where is a Lagrange multiplier enforcing the normalization constraint, and I have used the fact that is symmetric in to simplify the expression. Taking a functional derivative with respect to on both sides gives


where I have set the expression equal to zero in order to locate the extrema. Rearranging this equation, then taking a sum over on both sides gives


finally resulting in


We now have an equation involving two unknowns, and . We need a second equation in order to determine these two parameters. We may use the constraint equation to solve for , giving


Plugging this back in to our transfer matrix we find,


This tells us that is equivalent to the eigenvalue we sought, . As a check, we can compare Eqs. (52,81) in the unitary limit: , giving


for both Equations.

In Sec. II.2.5 we will discuss a simple formalism for determining the exact two particle spectrum in a box for any given scattering phase shift. This will allow us to eliminate certain finite volume systematic errors automatically. The transfer matrix method is also powerful because it gives us access to the entire two particle, finite-volume spectrum. When we discuss improvement in Sec. IV.2, we will add more operators and couplings to the interaction in order to match not only the ground state energy we desire, but higher eigenvalues as well. This will allow us to control the interaction between particles with non-zero relative momentum. To gain access to higher eigenvalues, the transfer matrix must be solved numerically, however, this may be accomplished quickly and easily for a finite volume system.

Iii Calculating observables

Perhaps the simplest observable to calculate using lattice (or any imaginary time) methods is the ground-state energy. While the two-body system may be solved exactly and used to set the couplings for two-body interactions, correlation functions for -body systems can then be used to make predictions. However, the transfer matrix for cannot in general be solved exactly, because the dimension of the matrix increases with particle number. For this reason we form instead -body correlation functions,




is a source for particles with spin/flavor indices , and a spatial wavefunction . For the moment the only requirement we will make of the wavefunction is that it has non-zero overlap with the ground-state wavefunction (i.e. it must have the correct quantum numbers for the state of interest).

Recall that a correlation function consists of insertions of the transfer matrix between source and sink. We can then expand the correlation function in a basis of eigenstates,


where is the overlap of wavefunction with the energy eigenstate , and is the th eigenvalue of the Hamiltonian. In the limit of large Euclidean time (zero temperature), the ground state dominates,


with higher excited states exponentially suppressed by , where is the energy splitting between the th state and the ground state. It should be noted that for a non-relativistic theory the rest masses of the particles do not contribute to these energies, so the ground state energy of a single particle at rest is , in contrast to lattice QCD formulations.

In this way, we can think of the transfer matrix as acting as a filter for the ground state, removing more excited state contamination with each application in time. A common method for determining the ground state energy from a correlation function is to construct the so-called effective mass function,


and look for a plateau at long times, whose value corresponds to the ground-state energy.

Once the ground state has been isolated, we can calculate matrix elements with the ground state as follows,