High Statistics Analysis using Anisotropic Clover Lattices: (II) Three-Baryon Systems

We present the results of an exploratory Lattice QCD calculation of three-baryon systems through a high-statistics study of one ensemble of anisotropic clover gauge-field configurations with a pion mass of m_\pi ~ 390 MeV. Because of the computational cost of the necessary contractions, we focus on correlation functions generated by interpolating-operators with the quantum numbers of the $\Xi^0\Xi^0 n$ system, one of the least demanding three baryon systems in terms of the number of contractions. We find that the ground state of this system has an energy of E_{\Xi^0\Xi^0n}= 3877.9\pm 6.9\pm 9.2\pm3.3 MeV corresponding to an energy-shift due to interactions of \delta E_{\Xi^0\Xi^0n}=E_{\Xi^0\Xi^0n}-2M_{\Xi^0} -M_n=4.6\pm 5.0\pm 7.9\pm 4.2 MeV. There are a significant number of time-slices in the three-baryon correlation function for which the signal-to-noise ratio is only slowly degrading with time. This is in contrast to the exponential degradation of the signal-to-noise ratio that is observed at larger times, and is due to the suppressed overlap of the source and sink interpolating-operators that are associated with the variance of the three-baryon correlation function onto the lightest eigenstates in the lattice volume (mesonic systems). As one of the motivations for this area of exploration is the calculation of the structure and reactions of light nuclei, we also present initial results for a system with the quantum numbers of the triton (pnn). This present work establishes a path to multi-baryon systems, and shows that Lattice QCD calculations of the properties and interactions of systems containing four and five baryons are now within sight.


I. INTRODUCTION
One of the ultimate goals of Lattice Quantum Chromodynamics (LQCD) is to calculate the properties and interactions of light-nuclei to high precision from first principles. While it is important to be able to explicitly demonstrate that nuclei emerge from QCD, the underlying motivation for this effort is to provide proof that LQCD provides a reliable theoretical tool with which to calculate highly complex low-energy strong interaction processes. With such a tool in hand, calculations of strong interaction systems for which experimental guidance is minimal, or absent, can be performed with confidence and with uncertainties that can be rigorously quantified. The interaction between three neutrons, which is an important input into many-body calculations of nuclei, provides an example of a quantity that is difficult to access experimentally, but which will be calculable to high precision with LQCD within the next decade.
The theoretical framework with which to determine the hadron-hadron scattering phaseshifts below the inelastic threshold from the volume dependence of the two-hadron energy levels in the lattice-volume was established a number of years ago by Lüscher [1,2]. This framework was used to extract nucleon-nucleon scattering lengths in quenched QCD (QQCD) [3] at unphysically large pion masses. Subsequent fully-dynamical LQCD calculations also used the Lüscher-method to extract nucleon-nucleon [4], and hyperon-nucleon [5] scattering lengths and phase-shifts 1 (from a single correlation function), albeit at unphysically large pion masses. The exponentially degrading signal-to-noise ratio encountered in the region of the correlation functions dominated by the ground-state, expected from the arguments presented by Lepage [13], severely limited the precision with which the scattering phase-shifts could be extracted from all of these calculations. The resources required to perform these calculations and the anticipated scaling of the statistical uncertainties in such calculations as a function of the pion mass [13] are sufficient to estimate the resources required to perform calculations of baryon-baryon (BB) scattering at the physical value of the pion mass to a given level of precision [10]. Given the smallness of nuclear physics energy-scales, typically a few MeV, compared to the chiral-symmetry breaking scale, Λ χ , or the pion mass, m π , it is guaranteed that a very large number of measurements will be required to achieve the necessary precision. At the physical pion mass, Ref. [10] estimates that ∼ 3 × 10 6 measurements will be required to determine the NN scattering length with a ∼ 50% uncertainty. Such estimates will be further refined as additional calculations at different pion masses, lattice-volumes and lattice-spacing are performed.
Recently, we have performed a high-statistics calculation of a number of single-hadron correlation functions [14] on an ensemble of the anisotropic gauge-field configurations generated by the Hadron Spectrum Collaboration [15,16] with a pion mass of m π ∼ 390 MeV, a spatial lattice spacing of b s ∼ 0.1227, an anisotropy ξ = b s /b t = 3.5 and a lattice volume of 20 3 × 128. The goal of the study was to "jump" an order of magnitude in the number of measurements performed to estimate correlation functions, and to explore the "new territory" that subsequently emerged. The baryon masses were extracted with precision at 1 Calculations of these same processes were subsequently performed in quenched QCD [6,7] and in QCD [8,9]. In these same works it was suggested that a phenomenologically useful baryon-baryon potential could be defined from LQCD calculations. However, flawed reasoning led to such a conclusion, in particular, the omission of the spatially-dependent two-body overlap factor which is present in the correlation functions [10,11,12]. the < ∼ 0.2%-level from the 292, 500 measurements performed on 1194 of these gauge-field configurations. A number of important and surprising observations were made in that work that have modified the path from LQCD to nuclei that we envisage. One of the most important aspects of that high-statistics work was that a detailed study of the signal-to-noise ratio in the single-baryon correlation functions became possible. The signal-to-noise ratio was found to be approximately independent of time for a significant number of time-slices prior to evolving toward the expected exponential degradation [13,14]. This window of "clean" time-slices is understood in terms of the the relative magnitude of the overlap of the single-baryon interpolating operator onto the single-baryon eigenstates, compared with the overlap of the corresponding interpolating-operator onto the lightest eigenstates (involving both meson and baryon-anti-baryon states) that contribute to the correlation function that governs the variance of the single-baryon correlation function. Given that the signal-to-noise ratio for a system containing more than one baryon is expected to scale (approximately) as the product of the signal-to-noise ratio's of the individual baryons (neglecting their interactions), this window of clean time-slices suggests that it may well be possible to calculate the energy-levels of systems containing a number of baryons in this lattice volume with these interpolating-operators.
In this work we present the first LQCD calculations of system comprised of three baryons. As the number of contractions required to form the correlation functions is naively N u ! N d ! N s !, one of the least computationally expensive systems 2 to explore is the one that couples to a source and sink of the form Ξ 0 Ξ 0 n. For simplicity, the product of the single baryon interpolating operators with the quantum numbers of the Ξ 0 and n are used in the calculations. This source will have non-zero overlap with states with strangeness four (s = 4), spin one-half (J = 1 2 ), and with a third component of isospin of one-half (I z = 1 2 ). Further, we present preliminary calculations of the pnn-system which will contain the triton if it is bound for this pion mass. We are presently unable to explore even the simplest system containing four baryons because of the computational resources required to perform the contractions (which are usually the least expensive component of a lattice calculation!), but see no reason why systems containing four and five baryons could not be explored in the near future.

II. LATTICE QCD CALCULATIONS
In this study, we employ a single ensemble of the n f = 2 + 1-flavor anisotropic clover gaugefield configurations that have been produced by the Hadron Spectrum Collaboration [15,16]. The technical details of the propagators computed on this ensemble are presented in Ref. [14] and we do not repeat them here. In the current calculation, the analysis is restricted to a slightly smaller data set, corresponding to an average on 218 randomly distributed measurements on each of 1191 configurations (a total of ∼260,000 measurements).
Each of the propagators calculated on the gauge-field configurations is used to determine two-point correlation functions, which for a single baryon have the form where H α (x, t) is an interpolating operator for the appropriate baryon state, e.g., for the proton H α (x, t) = abc u a,T C γ 5 d b u c,α (x, t) where C is the charge conjugation matrix. The Dirac matrix Γ is an arbitrary particle-spin-projector and the point x 0 is the propagator source point. The interpolating-operator at the source, H, is constructed from gaugeinvariantly-smeared quark field operators, while at the sink, the interpolating operator is constructed from either local quark field operators, or from the same smeared quark field operators used at the source, leading to two sets of correlation functions. For brevity, we refer to the two sets of correlation functions that result from these source and sink operators as smeared-point (SP) and smeared-smeared (SS) correlation functions, respectively. The correlation functions for the three-baryon systems have the form 3 , whereΓ is the tensor that projects onto the required angular momentum state. The same quark-propagators have been used in each baryon, and thus the source for each baryon is located at the same spatial point. More physically motivated sources and sinks involving spatial separations would likely improve the overlap onto the ground state in these systems, however this approach would be more computationally demanding and is not used in this exploratory work.
In the present work, we have restricted ourselves to the calculation of correlation functions for which each baryon is projected to zero-momentum at the sink, defining C H 1 H 2 H 3 ;Γ (t) = C H 1 H 2 H 3 ;Γ (0, 0, 0; t). Further, the optimal analysis of the three-baryon systems (with propagators from a single source) would have involved calculating the correlation functions associated with the different sink-smearing, SP and SS, for each baryon. Due to lack of computational resources we have restricted ourselves to the (SS) 3 and the (SP) 3 correlation functions, and have not calculated the "mixed" correlation functions, such as the (SS) 2 (SP) correlation function.
With two correlation functions associated with each set of quantum numbers in both the one-and three-baryon sectors, a linear combination of the pair of correlation functions can be constructed to produce a combined correlation function that more cleanly projects onto the lowest energy state in the lattice volume. One way to accomplish this is by hand, where one simply "looks" for the linear combination of correlation functions that has an effective mass plot (EMP) with the ground-state extending to the shortest time-slice. A refinement of this "brute-force" method is to use the matrix-Prony method presented in our previous paper [14], which we now review.
The two correlation functions, SS and SP, are sums of exponentials and satisfy the following matrix relation, where M and V are 2 × 2 matrices and y(t) is a column vector with two components corresponding to the two correlation functions. Eq. (3) implies then the correlation functions 3 A more complete calculation would generate correlation functions between sources and sinks that carry the same global quantum numbers, such as Ξ 0 Ξ 0 n → Ξ 0 Σ + Σ − , in order to identify all of the states in the lattice volume. For computational expediency we study only one combination of source and sink. are where q n and λ n = exp(m n ) are the eigenvectors and eigenvalues of the following generalized eigenvalue problem Given the two sets of correlation functions, the masses can be found by determining the matrices M and V that are needed in order for the signal to satisfy Eq. (3). Solving Eq. (5), leads to the eigenvalues and eigenvectors needed to reconstruct the amplitudes of each exponential in the correlation functions. It is straightforward to show that a solution for M and V is These inverses exist provided that the range, t W , is large enough so that the matrices in the brackets are of full rank. Once the eigenvalues, λ n and eigenvectors q n are determined, the amplitudes, A n , can be reconstructed using a fixed time-slice as a normalization point. The parameters t W and t J can be used to improve stability as investigated in Ref. [14]. The eigenvectors associated with the ground-state energy eigenvalue provide the linear combination of SS and SP correlation functions for which the plateau in the EMP sets in at the earliest time. As the eigenvectors can be determined at early time-slices, the degradation of the signal at later times seen in the eigenvalues of the matrix-Prony method, largely due to increasing fluctuations in the SS correlation function, is greatly reduced. This method is independently applied to the single baryon and the three-baryon pairs of correlation functions, to produce a single correlation function for each. For the present calculations, the relevant "diagonalized" correlation functions are where the coefficients η which is independent of time when the diagonalized correlation function is dominated by a single exponential. The effective mass plot (EMP) associated with the diagonalized nucleon correlation function is shown in fig. 1, and that associated with the diagonalized Ξ 0 correlation function is shown in fig. 2. Extended and clean plateaus are observed for both the nucleon and the Ξ 0 , as discussed in detail in Ref. [14].
The source and sink used to produce the Ξ 0 Ξ 0 n state is the product of interpolating operators that have good overlap onto the lowest-lying octet baryons, Ξ 0 and n. However, since the three-baryon eigenstates of the QCD Hamiltonian in the lattice volume, or in nature, are not simple products of single baryon eigenstates, this source and sink will couple (at some level) to all states with the corresponding quantum numbers. This will be the case, no matter how well the single hadron interpolating operators project to their respective ground states. We expect our correlation functions to have significant contributions from nearby states, such as thresholds (neglecting interactions) of 0.6893, 0.6893, 0.6858, 0.7008, 0.7008, 0.7008, 0.6933 in lattice-units, respectively. It is clear that the energy eigenstates in the lattice volume will be mixtures of the different states and from the above considerations, we expect to find four relatively close energy-levels. 4 In order to cleanly see these nearby states, multiple correlation functions formed from different sources and sinks, and more sophisticated analysis techniques will be required. As the goal of this work is not to provide detailed spectroscopy of such states, but to demonstrate the feasibility of studying such systems, we do not make efforts to identify the total isospin of the ground state, and are content with identifying what appears to be (with the current statistics) a single state. For simplicity, we refer to this state as Ξ 0 Ξ 0 n. The EMP's associated with the SS, SP and diagonalized Ξ 0 Ξ 0 n The uncertainties in this result correspond to a statistical uncertainty, a fitting uncertainty from the analysis presented above, and an additional fitting uncertainty from comparison to alternate analysis techniques using either multiple exponential fits or other Prony methods (see Ref. [14] for details). An additional uncertainty associated with the determination of lattice scale b s = 0.1227 (8) is not included. Fig. 3 also shows a plateau for a backward propagating negative parity state, which is consistent with a four-body state with the quantum numbers of Ξ 0 Ξ 0 nπ (which ultimately will allow for the calculation of pion interactions with multi-baryon systems). As the diagonalized correlation functions are dominated by their respective ground-states even at relatively short times, the energy-splitting between the ground-state of the Ξ 0 Ξ 0 n system, and that of the two Ξ 0 's and a neutron can be found efficiently by forming the ratio of the diagonalized correlation functions, C H (t), which at large times (for gauge-field configurations that are infinitely long in the timedirection) tends to an exponential that depends upon the energy-splitting δE Ξ 0 Ξ 0 n =  Fig. 4 shows the effective mass corresponding to G Ξ 0 Ξ 0 n (t), along with the correlated fit to the plateau region. The energy-splitting δE Ξ 0 Ξ 0 n is determined to be from fitting the time-interval from t = 21 to t = 35, using t J = 3, which is consistent with zero (the splitting is computed relative to the non interacting system of two Ξ 0 s and a neutron for convenience). It is very encouraging that the uncertainty in the energy-shift per baryon is ∼ 3 MeV, which is smaller than the binding-energy per nucleon in typical nuclei, B ∼ 8 MeV, and not significantly larger than the binding-energy per nucleon in the deuteron or triton at the physical values of the light-quark masses. The single energy-level fit to the EMP in fig. 4 has a χ 2 /d.o.f. = 2.0, indicating that there maybe additional structure in the correlation function. Including a second energy-level shifted by ∆E ∼ −0.004 lattice units might provide a better description of the EMP, and this would be consistent with the lower-energy state, Ξ 0 ΛΛ, that is expected to contribute to the four low-lying eigenstates in the lattice-volume. However, enhanced statistics are required to determine if this is, in fact, the case. At present, unlike the situation in multi-meson systems [17,18,19], the analytical tools are not in place to use the above energy shift and those of the associated two baryon systems to extract the parameters describing the relevant two-and three-body interactions. While the volume dependence of the simplest three-fermion systems has been studied in Ref. [20], the mixing we expect between four closely spaced states complicates the situation.

III. SIGNAL-TO-NOISE RATIOS
The Ξ 0 Ξ 0 n calculation is possible with our present resources because there are time-slices in the correlation functions for which the signal-to-noise ratio is approximately independent of time. This is a region of time-slices for which the correlation function that dictates the variance of the signal is not yet dominated by its "ground-state", which for the single-nucleon correlation function is three pions. Given the importance of this observation in Ref. [14], it is worth re-stating and expanding upon it here.
As argued by Lepage [13], correlation functions involving one or more baryons exhibit statistical noise that increases exponentially with Euclidean time. In the case of a single positive parity nucleon, the correlation function has the form where N α (x, t) is an interpolating field that has non-vanishing overlap with the nucleon, Γ + is a positive energy projector, and the angle brackets indicate statistical averaging over measurements on an ensemble of configurations. The variance of this correlation function is given by where all interaction energies have been neglected, and N is the number of (independent) measurements (distinct from the nucleon field operator N ). Therefore, at large times, the noise-to-signal ratio behaves as More generally, for a system of A nucleons, the noise-to-signal ratio behaves as at large times. The degradation of the signal-to-noise ratio on gauge-field configurations of finite temporal extent is exponentially more rapid than that given in Eq. (15) due to the presence of thermal states, as discussed in Ref. [14]. From the signals and variances that we have measured in the one-, two-and three-baryon sectors, it is clear that there is a suppression of the overlap onto the three-meson state from the N N source and sink (the variance correlation function of Eq. (13)), as encapsulated in the factor Z 3π . If Z 3π Z N N , Z 2 N there will be a number of time-slices, near the source of the correlation function, for which the noise-to-signal ratio behaves as and does not depend exponentially upon time, or the differences of hadron masses. The correlation functions we have constructed lead to an implicit suppression of Z 3π compared to Z N N and Z 2 N , due to the fact that the overlap onto the three-meson state, or any meson state, is strongly suppressed when the sinks N α (x, t) and N β (y, t) do not overlap within a volume approximately defined by the pion Compton wavelength. Therefore, summing independently over the volumes for N α (x, t) and N β (y, t) leads to a suppression factor that scales with the spatial lattice volume 5 as Z 3π /Z N N ∼ 1/Λ 3 L 3 . The pion mass dictates the weakest suppression, and hence we set Λ = m π for the estimates that follow (if the width, w, of the smearing of the source and sink is larger than the pion Compton wavelength, then Λ ∼ 1/w). Results consistent with this volume scaling have been found explicitly in calculations of single baryon energies using domain-wall fermions on MILC gauge configurations. By generalizing this argument to systems composed of A nucleons 6 where each interpolating field is projected to zero momentum, the noise correlation function is expected to behave parametrically as where we have made explicit the parametric dependence of the overlap factors on the baryon number and spatial volume. The dependence on A arises from the number of ways that N and N sink operators can overlap to form one or more three-pion contributions to the correlation function. Provided that the spatial volume is large compared to the Compton wavelength of the pion, m π L 1, there will be a range of time-slices in which only the last two terms in Eq. (17) are important. In this region, the signal-to-noise ratio in the multibaryon correlation function does not degrade exponentially faster that the signal-to-noise ratio in the single baryon correlation function. Instead, where Z and Z are O(1) ratios of overlap factors. Consequently, the signal-to-noise ratio starts degrading exponentially only after time-slice t noise , which has parametric dependence It is important to note that t noise depends only logarithmically on the number of baryons, and hence it is conceivable that plateaus may be found in the EMP's of systems containing four or more baryons with the current number of measurements if the contractions are performed. In order to investigate the signal-to-noise ratio in the correlation functions of interest, it is useful to form the effective noise-to-signal plot [14], in analogy with the EMPs. On each time slice, the quantity is formed, from which the energy governing the exponential behavior can be extracted via If the correlation function is dominated by a single state, and a single energy-scale determines the behavior of the noise-to-signal ratio, the quantity E S (t; t J ) will be independent of both t and t J . In fig. 5, the energy scales of the noise-to-signal ratio are shown for the nucleon and Ξ 0 . As discussed previously in Ref. [14], it is clear that for t 45 this scale is significantly greater than the asymptotic estimate of Lepage (the lowest horizontal line in each figure) because of thermal states involving propagation around the temporal extent of the lattice. It is also clear that even the simple Lepage scaling does not set in for many time-slices corresponding to a large window in which the signal is statistically clean. Figure 5 indicates that the suppression of mesonic intermediate states is stronger for the Ξ correlation function than for the nucleon correlation function, as evidenced by the energy-scale of the signal-to-noise ratio remaining small for longer times.
The energy-scale associated with the noise-to-signal ratio of the diagonalized Ξ 0 Ξ 0 n correlation function is shown in the left panel of fig. 6. While degrading exponentially, the signal-to-noise ratio of the three-body correlation function, is exponentially better in the The left panel shows the energy-scales associated with the signal-to-noise ratios for the Ξ 0 Ξ 0 n correlation function, as defined in eq. (21). The horizontal line corresponds to m N + 2m Ξ − 2m η − 5 2 m π , the asymptotic energy-scale in a lattice with infinite temporal extent. The right panel shows the difference between the signal-to-noise energy scales of the diagonalized Ξ 0 Ξ 0 n correlation function and that of the nucleon and twice that of the Ξ correlation function. plateau region t < ∼ 32, than expected based upon the arguments of Lepage, consistent with the expectations based upon the behavior of the signal-to-noise ratio of the nucleon and Ξ. The right panel of fig. 6 shows that the energy-scale associated with the signal-to-noise ratio in the Ξ 0 Ξ 0 n correlation function is consistent with the simple sum of the energy-scales from the single baryon correlation functions (within statistical uncertainties of the calculation).
In creating sources and sinks for correlation functions, a great deal of attention is paid to optimizing the overlap onto the states of interest. Variational techniques [22,23], the matrix-Prony method, and related approaches make use of sources and sinks with substantial, but different, overlaps onto the states of interest to enable a diagonalization to the eigenstates in the lattice-volume (up to exponentially suppressed contributions). For multi-baryon systems, the results of this work make clear that an equally important component of source and sink optimization is to minimize the overlap onto the lightest states that contribute to the variance of the correlation function. This will also be true for the extraction of the properties of excited single particle states. Figure 7 shows a comparison of the relative uncertainties (statistical and systematic uncertainties are added in quadrature and normalized by the mean value of the measurement) in the extraction of the ground state hadron energy for a selection of one-, two-and threebaryon systems using the measurements for which Ξ 0 Ξ 0 n contractions exist. We find that this quantity is approximately constant, due to there being a sufficiently large window of time-slices for which the signal-to-noise ratio does not degrade exponentially. The extent of this window is empirically seen to decrease with the number of baryons as shown in the lower panel of fig 7 where the matrix-Prony effective energies 7 of exemplary one-, two-and three-baryon systems are shown (detailed analysis of the two-baryon sector will appear in future work [24]). This result is consistent with the scaling anticipated in Eq. (19). 7 These EMPS and the central values of these extractions slightly differ from the extractions presented in figs. 1, 2, and 3, but are consistent. The standard EMP for the smeared-smeared pnn correlation function is shown (circles) along with the generalized EMP resulting from a 3-exponential Prony analysis [14] of the data (squares). The horizontal light band corresponds to three-times the nucleon mass.

IV. THE TRITON CHANNEL
One of the main motivations for the present work is to show that three-baryon (and beyond) calculations can be done at unphysical quark masses with present day resources. We have shown the results for channels that couple to Ξ 0 Ξ 0 n sources and sinks as their correlation functions are some of the least computationally expensive, requiring the calculation of only 288 Wick contractions. A significantly more computationally expensive, but physically more interesting, channel is that of the triton (pnn) for which there are 2880 Wick contractions. We do not have the computational resources available to perform the pnn contractions on all of the 260,000 measurements we have made, and to date have only performed ∼ 9, 200 measurements of the smeared-smeared correlation function. The generalized EMP resulting from a 3-exponential Prony analysis is shown in fig. 8, and given the relatively large uncertainties, we do not present a value for the ground-state energy. At the physical pion mass one expects to find a negatively shifted state corresponding to the triton. To conclude that such a signal corresponded to a bound state would require further studies showing exponentially suppressed sensitivity to the lattice volume in contrast to continuum states.

V. CONCLUSIONS
In this work, we have presented the first Lattice QCD calculations of a three-baryon state, focusing on a system with the quantum numbers of Ξ 0 Ξ 0 n. We find a ground state energy of E Ξ 0 Ξ 0 n = 3877.9 ± 6.9 ± 9.2 ± 3.3 MeV corresponding to an energy shift from the free three-baryon system of δE Ξ 0 Ξ 0 n = 4.6 ± 5.0 ± 7.9 ± 4.2 MeV. Our high-statistics analysis of the behavior of the signal-to-noise ratio of single-and multiple-baryon correlation functions indicates that there is a window of time-slices for which the signal-to-noise ratio does not degrade exponentially. This implies that multi-baryon correlation functions can be calculated in this time-interval with significantly less computational resources than previously estimated, and we demonstrate that this is indeed the case in the three-baryon sector. The signal-to-noise ratio does not depend exponentially upon the number of baryons in this timeinterval, however, the length of this window decreases logarithmically with the number of baryons.
We have focused only on the state(s) that couples to the Ξ 0 Ξ 0 n interpolating-operator simply due to limited computational resources. In the past it has been the case that gaugefield generation has required the majority of LQCD resources, but this is no longer true for precise calculations of baryonic observables. The resources required to perform the large number of measurements required for nuclear systems is significantly greater than that required for gauge-field generation. This situation will improve as more effort is put into algorithmic improvements for contractions, in the same way that the use of deflation [25] and other techniques have greatly reduced the resources required for propagator generation. Work in this direction is in progress. Given the observed behavior of the signal-to-noise ratio, we hope to be able to identify at least the ground state in systems of four and five baryons.
As the central goal for applications of lattice QCD to nuclear physics is the calculation of nuclei and their interactions, we have also presented the first calculations of the correlation function that would contain the triton if the calculations were at the physical pion mass. The statistics are very limited compared with the Ξ 0 Ξ 0 n correlation function, but it is encouraging to see that there is a clear plateau visible in the effective mass (within somewhat large uncertainties).
The increase by more than one order of magnitude in the number of measurements performed on a given ensemble of gauge-field configurations has given rise to a new understanding of how to pursue nuclear physics processes with Lattice QCD. Source and sink optimization involves two considerations (maximal overlap onto the baryon states and minimum overlap onto the mesonic states in the correlation function dictating the variance of the baryon correlation functions) to make optimal use of available resources. It is clear that, at unphysical values of the quark masses, high statistics calculations can be used to explore multi-nucleon systems (perhaps beyond A = 5) with present day resources.

VI. ACKNOWLEDGMENTS
We thank R. Edwards and B. Joo for help with the QDP++/Chroma programming environment [26] and D. B. Kaplan for discussions. KO would like to thank A. Stathopoulos useful discussions on numerical linear algebra issues and for his contribution in the development of the EigCG algorithm [25]. EigCG development was supported in part by NSF grant CCF-0728915. We also thank the Hadron Spectrum Collaboration for permitting us to use the anisotropic gauge-field configurations, and extending the particular ensemble used herein. We gratefully acknowledge the computational time provided by NERSC (Office of Science of the U.S. Department of Energy, No. DE-AC02-05CH11231), the Institute for Nuclear Theory, Centro Nacional de Supercomputación (Barcelona, Spain), Lawrence Livermore National Laboratory, and the National Science Foundation through Teragrid resources provided by the National Center for Supercomputing Applications and the Texas Advanced Computing Center. Computational support at Thomas Jefferson National Accelerator Facility and Fermi National Accelerator Laboratory was provided by the USQCD