High Statistics Analysis using Anisotropic Clover Lattices: (I) Single Hadron Correlation Functions

We present the results of high-statistics calculations of correlation functions generated with single-baryon interpolating operators on an ensemble of dynamical anisotropic gauge-field configurations generated by the Hadron Spectrum Collaboration using a tadpole-improved clover fermion action and Symanzik-improved gauge action. A total of 292,500 sets of measurements are made using 1194 gauge configurations of size 20^3 x 128 with an anisotropy parameter \xi= b_s/b_t = 3.5, a spatial lattice spacing of b_s=0.1227\pm 0.0008 fm, and pion mass of m_\pi ~ 390 MeV. Ground state baryon masses are extracted with fully quantified uncertainties that are at or below the ~0.2%-level in lattice units. The lowest-lying negative-parity states are also extracted albeit with a somewhat lower level of precision. In the case of the nucleon, this negative-parity state is above the N\pi threshold and, therefore, the isospin-1/2 \pi N s-wave scattering phase-shift can be extracted using Luescher's method. The disconnected contributions to this process are included indirectly in the gauge-field configurations and do not require additional calculations. The signal-to-noise ratio in the various correlation functions is explored and is found to degrade exponentially faster than naive expectations on many time-slices. This is due to backward propagating states arising from the anti-periodic boundary conditions imposed on the quark-propagators in the time-direction. We explore how best to distribute computational resources between configuration generation and propagator measurements in order to optimize the extraction of single baryon observables.


I. INTRODUCTION
One of the primary goals of lattice QCD (LQCD) is to calculate the properties and interactions of nucleons and, more generally, systems comprised of multiple hadrons. Precise exploration of the simplest multi-hadron systems has recently become possible with significant advances in computing resources, as well as through algorithmic and theoretical developments. The two-pion system π + π + is the simplest of such multi-hadron systems to calculate in LQCD, and current computational resources have allowed for a precise determination of the π + π + scattering length [1,2] at the ∼ 1% level. Recently, we have explored systems comprised of up to twelve π + 's [3,4] and also systems comprised of up to twelve K + 's [5] for the first time, allowing a determination of the three-π + and three-K + interactions. In general, a determination of the two-particle scattering amplitude, or multi-body interactions, with LQCD requires calculating the energy-eigenvalues of the system in the finite-volume [6,7,8,9]. The energy differences between the multi-particle energy-levels in the finite-volume and the sum of the particle masses determines the scattering amplitude or interaction. As processes of interest to low-energy nuclear physics are in the MeV energyregime, while the masses of the baryons and nuclei are in the GeV regime, the energy-levels in the volume must be determined to high precision to yield useful constraints and predictions for scattering amplitudes, phase-shifts and electroweak properties. Consequently, correlation functions of systems comprised of more than one hadron must be calculated with small statistical and systematic uncertainties ( 1 %) in order to provide useful information about low-energy nuclear interactions and nuclei.
The correlation functions associated with systems of baryons (and, more generally, states other than the pion) suffer from an exponential degradation of the signal-to-noise ratio as a function of time as argued by Lepage [10]. The scale that dictates this degradation is the difference between the total energy of the baryons in the system and half of the total energy of hadrons that contribute to the correlation function associated with the square of the interpolating operator for the baryon system. An example is provided by the two point nucleon correlator (N is an interpolating field with the quantum numbers of the nucleon) where, neglecting effects of the finite temporal extent which we discuss below (here Z and Z are overlap factors). Since M N − 3 2 M π > 0, the signal degrades exponentially in time with this exponent. Further, for multi-baryon systems, this exponent is scaled by the baryon number, B, and it consequently requires exponentially larger computational resources to calculate the properties of systems containing B > 1 baryons than a single baryon. In many regards, it is this signal-to-noise problem that distinguishes LQCD calculations of quantities typically of importance to nuclear physics from those typically of importance to particle physics.
The main motivation for our present work is to explore very high statistics calculations of the energy spectrum of B = 0, 1 correlation functions, quantifying the statistical scaling and identifying any issues that appear in the regime of precision calculations. More generally, we aim to assess the feasibility of extracting precise phase-shifts and multi-nucleon interactions from multi-baryon systems but we leave these discussions to subsequent work. Our focus is on the statistical scaling behavior of these measurements instead of on measuring physically relevant quantities. Consequently, we work with a single ensemble of gauge configurations that was produced by the Hadron Spectrum Collaboration [11] (the details of these configurations are discussed below). The analysis presented here enables us to identify a number of issues that will be important to LQCD calculations of quantities where exponentially degrading signal-to-noise ratios are a dominant concern: 1. While the classic argument of Lepage [10] concerning the behavior of the signal-to-noise ratios of baryon correlation functions seems to be on a solid theoretical footing, it has yet to be explored and verified through direct calculation. We examine the signal-tonoise ratios of the single hadron correlation functions in detail and present a modified version of the Lepage argument that incorporates the finite extent of the temporal direction of the gauge-field configuration, focusing on the case of quark propagators subject to anti-periodic temporal boundary conditions (BCs). Over large regions of the temporal extent of the lattice, the signal-to-noise ratio degrades exponentially faster than expected from the original Lepage argument, see Sec. VIII.
2. At present, and even more so in the past, the generation of gauge-field ensembles consumes most of the computational resources of LQCD calculations. 1 However, it is not clear what the optimal distribution of computational resources between gaugefield production and measurements (propagator calculations and contractions) is when one is interested in noisy quantities. To address this, we explore what can be accomplished by performing hundreds of measurements per configuration, and how precisely the baryon ground-state masses can be determined from an ensemble of 1194 configurations. We also study whether the measurements "saturate" after some critical number have been performed on one configuration (that is, exhibit little or no improvement in uncertainties after a critical number of measurements), finding for baryons that they do not, even up to ∼ 200 measurements per configuration. 3. Correlation functions that are determined to high precision are amenable to analysis with a variety of techniques, beyond those typically used successfully with low statistics data. On these anisotropic configurations, multiple (five or more) exponential fits to such correlation functions become stable as statistical fluctuations decrease, and the ground state energies can be extracted with high precision. We show that the generalized effective mass (EM) method, in which multiple energies are extracted from a linear system (a method developed by Gaspard Riche de Prony in 1795) also becomes useful for correlation functions with small uncertainties. As two (different but correlated) correlation functions are computed per species of hadron, this method is extended to construct the matrix-Prony method, which is found to be a very clean and effective tool for determining the ground-state energies. 4. While the correlation function generated by a single baryon interpolating operator will be dominated by the baryon ground-state at large times, it also contains contributions from all states that can couple to the operator. This includes multi-hadron states. The backward propagating component of the nucleon correlation function is dominated by the lowest energy negative-parity I = 1 2 state for the projection-operator we have applied to the correlation functions. By measuring the energy of this state, which is above the πN threshold and therefore is a continuum state, the phase-shift associated with the s-wave πN interaction is determined at this energy. The important point here is that this process contains disconnected diagrams, which are encoded in the gaugefield configurations, and do not require additional (of order the volume in number) calculations. 5. We also note that thermal states, while strongly suppressed, are seen in our high precision data. In these states, some part of the hadronic state propagates backward in time and can consequently manifest itself in the correlation function as an exponential with energy less than that of the zero temperature ground state. These contributions have amplitudes that are exponentially suppressed by the temporal extent of the configuration, but they can be extracted in certain temporal regions of the correlation function(s) where other components are also small. They can lead to pollution of the ground state signal.
The structure of this work is as follows. Section II introduces the details of the lattice calculations we perform, and in section III we discuss our expectations for the hadron spectrum on this ensemble. Section IV introduces the tools used in our analysis and presents detailed comparisons of the different methods we utilize. Following this, sections V, VI and VII present our main results for the pseudo-scalar mesons, ground-state baryons and negative-parity excited states, respectively. In sections VIII and IX we discuss the behaviour of noise in our measurements and investigate the scaling of uncertainties in hadron masses for varying numbers of gauge configurations and measurements. We conclude in section X. In subsequent works, we will address states with baryon number, B > 1.

A. Calculational Details
In this study, we employ an ensemble of the n f = 2 + 1-flavor anisotropic clover gauge-field configurations that are currently being produced by the Hadron Spectrum Collaboration [11]. These ensembles are being generated with dynamical tadpole-improved clover fermions and a Symanzik-improved gauge action (see Ref [11] for full details). All of the calculations that we present here were performed on a single ensemble of gauge-field configurations of size 20 3 × 128 with an anisotropy parameter of ξ = b s /b t = 3.5, a spatial lattice spacing of b = b s = 0.1227 ± 0.0008 fm, a pion mass of M π ∼ 390 MeV and a kaon mass of M K ∼ 546 MeV. The ensemble used in this study contains 1194 configurations taken at intervals of 10 trajectories, after allowing 1000 trajectories for thermalization. Ref. [11] provides a comprehensive analysis of autocorrelation times and thermalization. Some correlation is seen to be present between configurations separated by 30 trajectories.
The light and strange quark-propagators were computed using the same fermion action used in the gauge-field generation. We use the clover discretization of the fermion action as it requires significantly less computational resources than, for instance, the Domain-Wall discretization, in both the production of gauge-field configurations and in the calculation of quark-propagators, while retaining an O(b)-improved spectrum. Unlike the Domain-Wall discretization, the Clover discretization does not have a lattice chiral symmetry. At moderate lattice spacings, this may significantly impact the extraction of the properties and interactions of pions and kaons, but it is not expected to produce systematic uncertainties that are as significant in the properties and interactions of baryons(this remains to be verified and will not be addressed here). The quark propagators are calculated with anti-periodic BC's imposed on the timedirection and periodic BCs imposed on the spatial directions. As multiple propagators are calculated on each configuration, iterative solvers beyond the simple conjugate gradient algorithm can provide significant speed improvements. In particular, we employ the deflated conjugate gradient algorithm proposed in Ref. [12], and implemented in the Chroma lattice field theory library [13] as the EigCG inverter. In our typical production runs, we compute from 30 to 100 propagators in sequence, observing a factor of ∼ 7 improvement in inversion speed after the first few solves are used to deflate low-lying eigenvalues from subsequent inversions. Figure 1 shows the details of the propagators computed in this work. The histogram indicates the number of propagators computed on each of the 1194 configurations (averaged over four adjacent configurations for clarity) . The total number of propagators computed in this dataset is 292,500, an average of 245 propagators per configuration (we note that the maximum number of point-to-all propagators that could be computed on each of the configurations is ∼ 10 6 ) Each propagator is generated from a gauge-invariantly Gaussian-smeared source [14,15], on a stout-smeared gauge-field in order to optimize the overlap onto the ground-state hadrons. On each configuration, the locations of the propagator source points are chosen The separation between pairs of sources on a given configuration, defined to be the minimum distance between two sources, including the effect of the (anti-) periodic boundary conditions. The height of the bar at R = 0 corresponds to the total number of propagators.
randomly throughout the configuration. In fig. 2, we show a histogram of the 4d-separation, R, between the each pair of sources on each configuration. The shoulder at R ∼ 4 fm appears because of the (anti) periodic boundary conditions. The average source separation is R ∼ 2.9 fm and the source density is 3.43 fm −4 .

B. Correlation Functions
The propagators are used to compute two-point correlation functions which, for baryons, take 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. Similar correlation functions are used for the mesons. The interpolatingoperator at the source, H, is constructed from gauge-invariantly-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. We calculate the smeared-point and smeared-smeared correlation functions associated with the π + , K + (J π = 0 − ) mesons, and the N, Λ, Σ, Ξ (J π = 1 2 + ) baryons. For the baryons, the energy projectors Γ ± = 1 2 (1 ± γ 0 ) are used to project separately onto either the positiveor negative-energy (parity) states (one of which can be time-reversed and added to the other to improve statistics). Correlation functions associated with a given pair of interpolating fields are averaged over all sources on each configuration, producing one correlation function per interpolating operator pair per configuration.

C. Statistical Behavior
Before extracting results for observables, we analyze the statistical behavior of the measured correlators. As the computational cost of each measurement is much less than the computational cost of generating each configuration, performing multiple, O(10), measurements on each configuration is a practical way to cheaply reduce uncertainties and is an approach that has been used by many groups. Averaging the measurements on a given configuration produces a more accurate estimation of the correlation function on that configuration. A priori, one might argue that performing a significantly larger number, say O(100-1000), of measurements on a given configuration is an inefficient use of computing resources as the additional measurements will contain little or no new information and will not decrease the statistical uncertainty in the measurements of interest. This argument is likely true for configurations extending over small volumes. Physically, one expects there are a number of length scales associated with the possible "saturation" density of the measurements on a given configuration. As the lightest hadron is the pion, one expects the critical saturation density of measurements to depend parametrically upon the dimensionless quantity ρ/M 4 π , where ρ = N src /V where V is the four-volume, and N src is the number of measurements on the configuration. For a simple quantity such as the energy, E, of an eigenstate in the volume, one also expects to find a dependence upon ρ/E 4 . For instance, we expect a dependence upon the dimensionless quantity N src /(V M 4 N ) for the nucleon. Figure 3 shows the scaling of the uncertainty in the effective mass (the logarithm of the ratio of the correlator on adjacent time-slices) at one particular time-slice for the π + , K + , N and Ξ as a function of the number of measurements per configuration. This calculation was performed on 664 of the 1194 configurations in the ensemble, those for which we have made more than 200 measurements. The correlation functions, after being averaged over the sources on each configuration, were blocked in units of 10 configurations (100 trajectories), and the uncertainties in the effective mass (EM) on each time-slice were generated with the single omission Jackknife procedure. The π + and K + correlation functions clearly show deviations from statistical independence beyond ∼ 10 sources per configuration, and by 200 sources per configuration there is little to be gained by performing additional measurements on a configuration. In contrast, measurements of the baryon correlation functions are behaving as if they are statistically independent even with 200 sources per configuration. It is clear that the statistical uncertainties in the baryon correlators can be further reduced by performing even more measurements per configuration. These observations are consistent with the arguments regarding the critical source densities.
An alternative way to investigate this question is to consider the correlation between measurements of a correlation function from different sources on the same configuration. A natural quantity to consider is an extension of the standard autocorrelator to a source-tosource autocorrelator, χ src , defined by where C(t, c, s) is the correlation function of interest evaluated on time-slice t, configuration c and from source s and the function θ(s 1 , s 2 : R) is unity if the two sources are separated by a 4d-distance |s 1 − s 2 | < R. A nonzero value of χ src (R) indicates the presence of significant correlations over distances shorter than R. We have calculated χ src for a number of the correlation functions we analyze but find no sign of deviation from zero even for the case of the π + . This may in part be due to the poor statistics at small source-source separations (see Fig. 2). A further consideration is that for a given number of configurations, at some value of N src , the uncertainty in the measurements of a correlation function on a given configuration will become smaller than the uncertainty in the measurements over the entire ensemble. Once this limit is reached, it is pointless to perform further measurements without also increasing the ensemble size. Our measurements are far from this limit as is illustrated by fig. 4 where the uncertainties in the measurements of π + and N correlation functions on some individual configurations are shown as a function of the number of sources and compared to the overall uncertainty attained with the full ensemble.
An important consideration in generating high statistics measurements is the correlation between configurations. Ideally, enough trajectories separate each gauge-field in the ensemble so that they are statistically independent to the precision of the calculation of interest. The degree of correlation between configurations dictates the number of measurements that should be performed on a given set of configurations before it is more cost effective to enlarge the ensemble. In fig. 5 we show the uncertainty at given time-slices of the EM for the π + , K + , N and Ξ as a function of the number of gauge-fields on which 50 measurements are performed. The configurations are maximally separated in Monte-Carlo time, but an increasing number of configurations means a reduced separation in Monte-Carlo time between each configuration. The curves in fig. 5 correspond to what is expected for statistically independent configurations. The 100 maximally separated configurations are separated by 80 trajectories, the configurations separated by 20 trajectories appear to be contributing as one expects for statistically independent configurations (assuming that those separated by 80 trajectories are statistically independent). This is consistent with the hadronic autocorrelation times measured on sets of configurations similar to this ensemble in Ref. [11], τ π ∼τ N ∼ 40.

D. Computational Costs
These calculations required significant computing resources to perform; the total cost of the measurements was approximately 7 × 10 6 JLab 6n cluster node hours (this is an older machine with a dual core 3 GHz Pentium D processor per node) distributed over various computational facilities. To put this into context, the generation of the gauge-field configurations required approximately one-third of this time [16]. Our calculational method makes use of hadronic building blocks (partly contracted sets of propagators) which are extremely useful for contracting multi-particle correlation functions but inefficient for single hadron correlation functions; only 4 × 10 6 JLab 6n cluster node hours were directly relevant to the calculations presented herein. Nevertheless, it seems that in order to achieve the level of precision of the B = 1 correlation functions presented here, propagator generation rather than gauge-field generation is the most computationally intensive component of the LQCD calculation. However, even this will be superseded by the calculation of the contractions that are required for systems involving more than two baryons (the subject of future work).

III. EXPECTED SPECTRA
The form of the correlation functions that are expected to emerge from these calculations is a textbook discussion, but is now becoming more relevant as advances in the field are enabling more complicated processes to be explored, such as scattering, excited states and multi-hadron interactions. Additionally, the accurate statistical sampling we perform in the current work brings to light features that have been safely neglected in the past. A discussion of the impact of the boundary conditions (here we use anti-periodic temporal BCs for the quark fields) on multiple meson correlation functions that were used in a recent calculation of K + K + scattering can be found in Ref. [5], and a more detailed derivation for a two particle system can be found in Ref. [17].
Baryon correlation functions are somewhat different, as the interpolating operator for the single nucleon, for instance, can couple not only to the N, but also to a state containing the N and an even number of π's, to a p-wave ΛK state, to a p-wave ΣK state,to a p-wave N π state, and to any other state with the same quantum numbers as the nucleon. Further, it can also couple to backward propagating negative-parity states, such as an s-wave N π. Finally, the single nucleon interpolating operator can couple forward and backward propagating hadronic states (these are thermal states as they exist only because of the finite temporal extent (temperature) of the configuration), an example being a forward propagating N and a backward propagating π or vice versa. These states are simply illustrated by an example shown in fig. 7, a N π thermal state. Here the finite temporal extent of the configuration is indicated by the vertical lines (these should be (anti-)identified). The two grey regions correspond to the source and sink interpolating field. In the case depicted, the interpolating field at the source is N = (uCγ 5 d T )u and that at the sink is N = (u T Cγ 5 d)u, suppressing spin and color indices. For the usual zero temperature ground state, the source produces three valence quarks and the sink annihilates three valence quarks. In the thermal state depicted, the source (right grey region) produces two valence anti-quarks and a valence quark (solid lines) while also producing, via gluonic interactions, a sea quark-antiquark pair (dashed line). The three anti-quarks between the grey regions combine to form an antinucleon propagating as exp(−M N (T − t)) where t is the separation between the source and sink and we ignore excited states for simplicity. The quark-anti-quark pair propagating around the temporal boundary (since two quarks propagate, the boundary appears periodic at the hadronic level) contribute a factor of exp(−M π t) where T is the temporal extent of the configuration. The resulting contribution to the two point correlator is then corresponding to a state with energy M N − M π in the observed spectrum.
In the case of the correlation function resulting from a single nucleon interpolating opera-tor in the A 1 representation of the cubic group, one expects to see a state with energy equal to the N mass, M N , and also states (a subset of all states) with energy M N − 2M π − δE ππ , M N + δE N π , and M N + 2M π + δE N ππ , where δE ππ is the interaction energy of two π's in an I = 0 state, δE N π is the interaction energy of the N π system, while δE N ππ is the interaction energy of the N ππ system. Particularly disturbing is the state with energy M N + δE N π corresponding to N π moving forward in time and a π moving backwards in time, that conspire to produce a state with an energy that differs from the nucleon mass only by the N π interaction energy. Such states will be exponentially suppressed by the temporal extent of the configuration, however, accurately disentangling such states from the zero temperature ground state will ultimately require calculations on ensembles of gauge-field configurations with different temporal extents.
It is important to realize that thermal states are not simply a curiosity that can be safely ignored. As we shall see in Section VIII, they dominate the statistical uncertainty of baryon correlation functions at large times, providing deviations from the naive form of the signal-to-noise ratio. The amplitudes of these states are exponentially suppressed by the temporal extent of the configuration times the mass of the backward going hadronic state. Consequently, the most important thermal states involve backward propagating pions, and, to suppress these states, the product M π T must be large. As the chiral limit is approached this will become more and more difficult since, in the limit, it is impossible to separate any particular state from itself and any number of pions.

A. Multi-Exponential Fits
The high statistics accumulated for this work allows us to perform stable multiexponential fits using a standard χ 2 minimization. In this section, we explore the determination of the ground and excited states as a function of several variables; the number of exponentials used in the fit function, N exp , the range of the fit, R, the number of sources per configuration, the number of configurations, the blocking time τ block , and the (effective) anisotropy, ξ eff . We present details of our fits for the Ξ, using a correlated fit to the smeared-smeared and smeared-point correlation functions.
To begin, we performed combined multi-exponential fits by minimizing where y s (t) are the lattice measured correlation functions, s = [SS, SP ], and Cov is the covariance matrix between both time-slices and correlation functions. The fitting functions used are, where C SS (C SP ) denotes the smeared-smeared (smeared-point) correlation function. To perform these fits we start with a single exponential and perform the correlated fit to Z S 0 , Z P 0 and E 0 . A selected set of best fit parameters from this fit are used as initial estimates for the two exponential fits. This is performed recursively by taking the best fit results from the N exponential fit as an initial estimates to the N+1 exponential fit. With this strategy, successful minimizations with up to six exponentials have been performed. However, with the inclusion of the fifth and higher exponentials, the minimizer performs poorly, and often returns two masses that are degenerate within their uncertainties. Furthermore, as discussed in detail in the previous Section, the expected spectrum of states on these anisotropic configurations is such that the resulting masses for the excited states are likely averages of nearby energy levels, see also Section IV E below for demonstrations of this. For these reasons, we are only confident in the ground state energies extracted in these fits. However, the number of exponentials used in a successful minimization plays an important role in minimizing the fitting systematic uncertainty.
The extracted mass of the Ξ as a function of the number of exponentials in the fit form is detailed in Table I. With the high statistics in this study, fitting a single exponential yields a statistical uncertainty of less than 0.2%, with a slightly larger fitting systematic uncertainty, however, 50 time-slices must be discarded because of excited state contamination. For our multi-exponential fits, the fitting systematic uncertainty is defined to be the standard deviation of all successful fits in a given minimization. To define a successful fit, we take a fixed length in time, R ≡ t max − t min , and a fixed set of initial parameters, and keep all fits with an integrated probability distribution Q > 0.1 while varying t min . 3 One observes that the statistical and systematic uncertainties are not further reduced by including more than three exponentials in the fit. The resulting ground state mass of the Ξ as a function of the t min used in the fit is shown in upper panel of fig. 8 (in a style similar to an effective mass plot) with the color and symbol shapes indicating the number of exponentials in the fit. The extraction of the nucleon mass is also shown in the lower panel. Increasing the number of exponentials in the fit, N exp , allows the t min -interval over which the ground-state energy is seen to plateau to be brought closer to the source where statistical uncertainties are much reduced.
The full set of measurements have been used to generate the fits presented in Table I, and the correlation functions from configurations nearby in Monte-Carlo time have not been TABLE I: Multi-exponential fits to smeared-smeared and smeared-point Ξ-correlation functions as a function of the number of exponentials, N exp . The number of successful fits, N suc has been defined to be fits with Q > 0.1 for a fixed time-length R with a fixed set of initial parameter values. The time window listed is taken as a representative example of a good fit. The first uncertainty is statistical while the second is taken from the standard deviation of successful fits, as defined above. blocked. Blocking is known to be important, since for correlated configurations, unblocked correlation functions can lead to underestimates of the true uncertainty. For hadronic quantities, we expect that the ensemble we have used has an auto-correlation time of about 40 Monte-Carlo time-steps [11]. Our calculations have been performed on configurations separated by only τ = 10. Several different fits were performed to determine the effects of blocking on our multi-exponential fits, the results of which are collected in Table II. To normalize these fits, a fitting interval, with a range of R = 30, was determined from the N exp = 3 exponential fits. For these fits, the blocking time τ block has no detectable impact on either the statistical or systematic uncertainties (τ block = 1 corresponds to no blocking, while τ block = 10 corresponds to blocking every 10 configurations). The range of time used in the fits also plays an important role in minimizing the uncertainty. The resulting fits for short and long time-ranges used in three and four exponential fits are shown in Table III. While the range does not have a significant impact on the statistical uncertainty, it does significantly reduce the systematic uncertainty in the fit. To have such long ranges of statistically useful time-slices, the anisotropy ξ = b s /b t , which is 3.5 for this ensemble, is crucial. We have not performed calculations with a different anisotropy (including isotropy), but this can be qualitatively studied by constructing correlation functions using only every second or every third time slice, with an effective anisotropy of ξ eff = 1.75 and 1.17, respectively. In the lower section of Table III, we display fits of 1, 2 and 3 exponentials to these reduced sets of measurements. This reduced anisotropy has a significant impact on the resulting uncertainties, particularly for ξ eff = 1.17. We were unable to find successful four exponential fits, and the number of successful fits with 1, 2 and 3 exponentials has been reduced. Furthermore, the smaller number of time-slices in the same physical extent, reduces our ability to control the systematics of the fits.
Finally, the impact of the number of sources, N src and number of configurations N cfg , on the uncertainties in the extracted mass of the Ξ has been explored, the results of which are collected in Table IV. With N src = 50 or N cfg = 597, the statistical uncertainties with N exp = 3 exponential fits are the same as with the full set of measurements. 4 However, in both cases the systematic uncertainty is larger than that of the full set of measurements. The corresponding dependence of more complicated multi-particle observables on the number of configurations and sources are under investigation.
For multi-exponential fits, it appears that the most important feature in controlling the uncertainty the ground state is the number of exponentials with which a successful min- imization can be performed. Neither the statistical nor systematic uncertainties improve beyond the inclusion of three exponentials in the fits. To have confidence in the N exp = 3 exponential fits, the anisotropy is found to be essential. A quantitative exploration of the effects of the anisotropy on the stability of multi-exponential fits is desirable, but this would be a very costly numerical endeavor. With three or more exponentials, the fitting range,

B. Generalized Effective Mass Plots
Correlation functions on an ensemble of configurations of infinite extent in the time-direction become dominated by a single exponential at large times with an argument that is the energy of the ground state of the system, It is conventional to define the effective mass (EM) from the logarithm of the ratio of the correlation function on adjacent time-slices. It is also possible 5 to form a more general EM from time-slices separated by t J > 1 For exponentially decreasing signals with time-independent noise, this will naturally reduce the statistical uncertainty in the EM and improve the extraction of energy-eigenvalues as it increases the "lever-arm" of the exponential. In such a case, the uncertainty in M eff (t J ) in Eq. (9) will decrease as 1/t J . Simple correlation functions involving π's have timeindependent uncertainties, but this is not the case for baryonic correlation functions, whose relative uncertainties grow exponentially with time. We explore the improvements to baryon EMs, and ultimately the extraction of baryon masses and the energy-eigenvalues in the volume, that result from t J > 1. In fitting an energy to an EM (and other generalizations), either the Bootstrap or Jackknife procedures are used to generate the covariance matrix associated with the time-slices in the range of the fit. 6 This covariance matrix is then used to form the χ 2 /dof, which is minimized to determine the energy, and then explored to determine the uncertainty in this energy. The statistical uncertainty is obtained by finding values of the fit parameters where the χ 2 function attains a value of χ 2 min + 1. To demonstrate the impact of t J > 1 for baryon EMs, we examine the smeared-smeared correlation function of the Ξ-baryon. Figure 9 shows the EMs obtained with t J = 1 (left panel) and t J = 10 (right panel). The scatter of the effective mass from time-slice to time-slice is significantly reduced with t J = 10 compared with t J = 1, allowing for a clear identification of the time range over which it is reasonable to extract the (ground-state) mass of the Ξ. Therefore, the systematic uncertainty associated with the fitting range in the EM is reduced. The statistical uncertainty in the mass of the Ξ extracted from the EMs with the two different values of t J , when fit over the time time-slice interval, are however very similar, as can be seen in the resulting fits to time-slices t = 48 to t = 58, The first uncertainty corresponds to the statistical uncertainty in the mass determined from the χ 2 /dof minimization, while the second corresponds to systematic uncertainty associated with the fitting interval. The systematic uncertainty of this fit is determined by varying the fitting interval at each end by 0, ±1, ±2 time-slices, performing a χ 2 /dof minimization over each interval and taking half of the spread of the extracted masses. Alternative procedures such as using fits to rolling windows of time-slices within the fitting interval return similar uncertainties.
It is interesting to explore how different values of t J modify the form of the covariance matrix that is input into the χ 2 /dof minimization. The covariance matrices associated with the time-interval t = 48 to t = 51 from these two EMs are shown in Eq. (11). They are quite 20 different, with the distant off-diagonal elements becoming more significant for increasing t J .
In this comparison, it is important to note that the two extractions make use of different parts of the correlation function. The t J = 1 fit uses five time-slices, while the t J = 10 fit uses eight well-separated time-slices. The EMs from the smeared-point Ξ correlation functions with t J = 1 and t J = 10 are shown in fig. 10. The scatter in the EM for t J = 1 is substantially less than for the smeared-smeared correlation function, as the overlap of the interpolating operator onto the ground-state is larger. However, the overlap onto excited states is even larger, and the ground state component of the correlation function does not become dominant until later times, increasing the fitting systematic uncertainty in the extraction of the ground-state mass.

C. Prony's Method/Linear Prediction
As the signal-to-noise ratio of baryon correlation functions degrades exponentially with time, it is important to extract the ground-state signal (or excited state signal if that is the state of interest) from a range of time-slices starting at the earliest possible time. Significant effort has been placed into determining interpolating operators that maximize the overlap onto the ground-states of the baryons in order to facilitate this. Further, there has been significant effort put into using the variational method [20,21], for which the correlation functions resulting from a number of hadronic interpolating operators are diagonalized on each timeslice to give the eigen-energies with the appropriate quantum numbers. A few years ago, Fleming suggested that generalizing the EM method to two or more exponential functions might be useful in LQCD analysis based on findings of NMR spectroscopists [22] 7 . At that time, we explored this technique with sets of correlation functions that were available to us at that time, and found the method was quite unstable to the statistical fluctuations in those measurements. More recently, Lin and Cohen [25] compared this method favorably to the variational approach. Given the small statistical uncertainties in the correlation functions we are presently considering, and the reduction in the systematic uncertainties achieved with t J > 1, we return to explore this technique.
In LQCD, two-point correlation functions have the form where t denotes the time-slice (time-slices are implicitly taken to be evenly-spaced). It follows from Eq. (12) that where the integer n is the generalization of t J to the case of a multi-exponential function.
In order to determine the k coefficients C i , k equations are required to be formed from the measured correlation function. Given the C i , the roots of and in particular the α's, provide the energies of the states contributing to the correlation function.

One Exponential : The Standard Effective Mass
In the case of k = 1, where the correlation function is assumed to be a single exponential, and taking n = 1, and the usual expression for the EM follows trivially.

Two Exponentials
In the case of two exponentials in the correlation function, the most general pair of equations that can be used to extract the two effective masses is where q 1 is an arbitrary integer off-set between the two equations. Inserting the values of the calculated correlation function allows for an extraction of C 0,1 on each time-slice, j. These coefficients are then inserted into to recover a numerical value for e −nα . By choosing n = m = 1, the expressions of Fleming [22] are recovered. In order to optimize the two-exponential extraction, a search over values of the pair (n, q 1 ) must be performed. A further systematic uncertainty can be assigned from this choice. The ground-state extracted from the smeared-smeared Ξ correlation function with n = q 1 = 5 is shown in fig. 11. It is clear that the ground-state signal can be isolated from the correlation function for a large number of time-slices, many more than using the single exponential EM ( fig. 9) alone. We have shown the fit to the ground-state result between time slices t = 30 and t = 60. The lower-limit of the time interval was chosen to be within an interval for which χ 2 /dof < 1. Extending the fit interval to lower time-slices gradually increases the χ 2 /dof, as shown in the right panel of fig. 11, indicating contamination from higher energy states. The upper-limit of the fitting interval was chosen to be in the region for which backward propagating states (due to the anti-periodic BC's in the time-direction) were not visible in the EM (or in the χ 2 /dof). The ground-state Ξ mass we extract from this 2-exponential analysis is where the first uncertainty is statistical and the second is the fitting systematic (as defined previously). The statistical uncertainty in the 2-exponential extraction is significantly smaller than that obtained from the one-exponential analysis (Eq. (10)). This is due to the substantially increased number of time-slices in the ground-state plateau in the generalized EM.
One aspect of this method that is less appealing is the ambiguity in the association of the two roots that result from Eq. (17) to the two states on different time-slices and on different jackknife/bootstrap ensembles. This (mis-)identification issue is the cause of the anomalously large uncertainties at time-slices 31,. . . ,36 in fig. 11 -this should not be interpreted as variance of the signal for the ground state. Additionally, on different timeslices and jackknife/bootstrap ensembles, this method can, and likely will, select different terms in Eq. (12) particularly for the sub-dominant excited state, adding additional artificial variance to the signals for particular energy eigenstates. Consequently, the extracted second state is not physically meaningful.

Three and More Exponentials
The generalization of the method to arbitrary numbers of exponential functions is straightforward. In the case of three exponentials, inserting the values of the calculated correlation functions, with q 1 = q 2 = 0 allows for an extraction of C 0,1,2 on each time-slice, t. Again, these coefficients can be extracted uniquely in terms of the G(t) due to the fact that the system is linear. These coefficients C 0,1,2 are then inserted into to recover a numerical value of e −nα . Analysis of a given correlation function involves searching for the values of the triplet (n, q 1 , q 2 ) that optimizes the extraction (in each equality in Eq. (19), different n can be used). In this case, statistical fluctuations occasionally result in complex roots of Eq. (20) on a particular Jackknife or Bootstrap ensemble. At present, we simply omit these contributions in our analysis. Such complex roots correspond to an oscillatory solution and arise from short distance noise in the correlation function (or nearly degenerate states in the spectrum), and are a well-known issue with the simple Prony method. More advanced methods [26] can mitigate this issue, but do not result in improved extractions of the ground-state so we do not discuss them in detail. The ground-state energy extracted from the smeared-smeared Ξ correlation function with n = 10, q 1 = 3, q 2 = 6 is shown in fig. 12. It is clear that the ground-state signal is extractable The left panel shows the ground-state extracted from the smeared-smeared Ξ correlation function with a 3-exponential Prony determination with n = 10, q 1 = 3, q 2 = 6, and the correlated fit to the time-slices between t = 10 and t = 52. The inner (yellow) region corresponds to the statistical uncertainty, while the outer (red) region corresponds to the statistical and fitting systematic uncertainties combined in quadrature. The right panel is the χ 2 /dof for fits between the time-slices t = t min (the horizontal axis) and t = 52.
from time-slices even closer to the source than with two-exponential analysis. In fig. 12, the fit to the ground-state between time slices t = 10 and t = 52 is shown. The extracted mass is M Ξ = 0.24124 ± 0.00032 ± 0.00034 , with the statistical uncertainty being slightly less than in the two-exponential analysis. It is important to realize that this level of precision corresponds to a statistical uncertainty of ∼ 2 MeV in the Ξ mass. We have successfully applied the four-and five-state Prony method to our data but no improvement is seen beyond the three-exponential extractions.

D. Multi-Correlation Function Prony Method
There are a number of extensions of the Prony method that exist in the literature (see for example [26]), some of which we have investigated in detail. For the correlation functions we have in hand, these extensions do not significantly improve on the standard Prony method. Typically, these methods are applied in cases where only a single set of measurements is available. However, we have two sets of correlation functions (smeared-smeared and smearedpoint) whose energy spectra are identical in the limit of a large number of configurations. It is straightforward to generalize Prony's method to include both correlation functionsthe matrix-Prony method. This form leads to a further reduction in the uncertainty of the extraction of the energy eigenvalues. A similar approach, has been briefly discussed in Ref. [27].
Assume we have N (N = 2 in our case) correlation functions from which we want to extract the energy levels. If these correlation functions are a sum of exponentials they satisfy the following recursion relation, where M and V are N × N matrices and and y(t) is a column vector of N components corresponding to the N correlation functions. Eq. (22) implies then the correlation functions are where q n and λ n = exp(m n t J ) the eigenvectors and eigenvalues of the following generalized eigenvalue problem M q = λV q .
Given the N 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. (22). Solving Eq. (24), then leads to the eigenvalues λ n = exp(m n t J ) and the eigenvectors q n needed to reconstruct the amplitudes with which each exponential enters the correlation functions. A simple solution can be constructed as follows. First note that Clearly, a solution for M and V is where these inverses exist provided that the range, t W , is large enough to make the matrices in the brackets full rank (t W ≥ N −1). In our case with two exponentials the range has to be two for achieving full rank. Once the eigenvalues, λ n and eigenvectors q n are determined, the amplitudes, C n , can be reconstructed using t as a normalization point. The shift parameter t J can be used it improve stability. The above solution is equivalent to determining M and V by requiring that is minimized.
To go beyond extracting two states, one can construct and solve a second order recursion relation. The minimization condition of Eq. (27), augmented to contain the second order terms in the recursion, can be used to determine the unknown matrices. The resulting eigenvalue problem is a second order nonlinear generalized eigenvalue problem which is straightforward to solve. However, to isolate the ground-state, which is our present focus, the two state model is sufficient and we do not pursue this further.
To demonstrate how this method works, we return to the Ξ correlation function discussed above. Figure 13 shows the generalized EMP for the Ξ mass as a function of time determined with a N = 2 matrix-Prony extraction, using both the smeared-smeared and smeared-point correlation functions. The inset shows the second extracted state in addition to the ground state. The extracted value of the Ξ mass, determined by fitting in the time interval t = 11 to t = 50, is M Ξ = 0.24097 ± 0.00025 ± 0.00003 , χ 2 /dof = 0.81 .
The EM of the dominant state in fig. 13 plateaus around time-slice t = 10, and is welldefined over a large interval. In addition to being somewhat more visually appealing than the previous Prony analyses of single correlation functions, this method provides the smallest uncertainties, particularly for the fitting systematic. In our final extractions of baryon masses, our EM analysis will use the matrix-Prony method. This method yields ground state energies that are in complete agreement with those from the other methods discussed. The generalized EMs from the matrix-Prony method are consistently clean, and the quality of fits are uniformly good for the ground state. Since they involve only one fit parameter, one can easily assess the quality of the fits. The procedure for fitting parameters and determining their statistical uncertainty has been described in Section IV B. Systematic uncertainties are calculated by performing fits over rolling windows of time-slices within the quoted overall range and looking at the standard deviation of the central values of those fits. This is combined in quadrature with a further systematic uncertainty that is generated by sampling a large range of possible values of t J and t W and taking the standard deviation of the central values of the resulting fits. The generalized EMP for the Ξ extracted with the matrix-Prony method for a variety of values of t W and t J can be seen in fig. 14.

E. Prony-Histograms
In extracting the energies of the states through the Prony procedure, a set of roots are produced on each time-slice for each member of the Bootstrap or Jackknife ensemble. In general these roots are real and there is an ambiguity in associating the roots with energylevels in the finite volume (only the single particle masses are approximately known). In order to aid identification of energy-levels it is useful to form histograms of the complete set of roots generated through the Bootstrap procedure. The simplest histogram is formed by accumulating all of the roots obtained on a sub-set of, or all of, the time-slices over all bootstrap/jackknife ensembles. The dominant components of the correlation function will appear as well-defined peaks in the histogram. In most cases, this histogramming procedure produces very similar results when either two, three or four exponential Prony, or matrix-Prony analyses are used. Only atypically do the higher exponential analyses reveal a clean state that is not present in the two-exponential analysis. Additionally, since our baryon correlation functions are asymmetric in time because of the parity projectors used in Eq. (2), noise is reduced in these histograms by separately accumulating the roots over the two half configuration. As expected, the excited states have a larger presence in the smeared-point correlation function. An example histogram is shown in fig. 15, corresponding to the Ξ correlation function analyzed in fig. 13. There is one clear peak in the histogram, corresponding to the Ξ ground-state and one broad structure at higher mass, which the histogram suggests is likely to be a collection of closely spaced states that currently are not resolvable. This interpretation of the excited state is consistent with expectations for the Ξ spectrum and with the instability of the extractions of the first excited state in the exponential fits discussed earlier.

V. MESON SPECTRUM
The π + correlation functions do not suffer from exponential signal-to-noise degradation for configurations of infinite temporal extent (in Section VIII, we will find that this is not true for a finite time-direction even for the pion). As a result, they can be calculated with small statistical uncertainty on each time-slice, as shown in fig. 16. As the EM for the π + does not exhibit a plateau, the π + mass is determined by fitting cosh M π (t − T 2 ) to a (large) number of time-slices of the correlation function. Performing a double cosh fit to time-slices t = 21 to t = 41 yields M π = 0.06936 ± 0.00012 ± 0.00005 , where the first uncertainty is statistical and the second is fitting systematic. The statistical uncertainty in the mass is determined with the Jackknife procedure, and the fitting systematic is determined by varying the fitting interval over a reasonable range. The second set of peaks that are visible in the histogram are at an energy consistent with the I = 1 KKπ state (with a threshold at M π + 2M K ∼ 0.2636) that can couple to the source that produces a single π + . With even greater statistics, the energy of this state could be calculated with enough precision to extract the I = 1 2 Kπ and I = 0 KK scattering lengths and the I = 1 KKπ three-body interaction. An expression for the energy-levels of this system in a finite volume in terms of the KK and Kπ scattering amplitudes and various three-body interactions has recently been derived [28] and would be useful in analyzing this state. There is no clear peak that can be associated with I = 1 πππ, which one would naively thought would have been present. It must be the case that the source does not couple with any appreciable strength to this state. The EM associated with the smeared-point kaon correlation function is shown in fig. 17, along with the bootstrap-Prony histogram. Despite the appearance of the EM, no plateau is found in the EM, and the kaon mass is extracted by fitting cosh M K (t − T 2 ) to a number of time-slices of the correlation function. Performing a double cosh fit over the time-slices between t = 29 and t = 49, yields a K + mass of M K = 0.097016 ± 0.000099 ± 0.000033 , χ 2 /dof = 1.01 .
The excited state(s) that are seen in the histogram in fig. 17

VI. GROUND-STATE BARYON SPECTRUM
With the methodology we have presented in Section IV, we are in a position to extract the masses of the lowest-lying octet baryons. The Ξ correlation functions have been used extensively to demonstrate the strengths and weaknesses of the various methods, with the resulting mass extraction given in Eq. (28), and we do not repeat that discussion here. The matrix-Prony method applied to the smeared-smeared and smeared-point correlation functions associated with the Σ, Λ and N produces the Prony-histograms and generalized EMs shown in fig. 18, fig. 19, and fig. 20. Fitting the Σ EM between time-slices t = 12 to t = 47 yields Σ mass, Fitting the Λ EM between time-slices t = 12 to t = 52 yields Λ mass, Finally, fitting N the EM between time-slices t = 11 to t = 40 yields N mass, The results of the best extractions of the ground-state baryon masses using multiexponential fitting and the matrix-Prony method, which give consistent results for each species of baryon, are collected in Table V. These results are completely consistent within their uncertainties, giving us confidence that our extractions are correct.

VII. NEGATIVE-PARITY EXCITED BARYON STATES
The interpolating operators that produce even-parity baryons moving forward in time also produce negative-parity partners moving backwards in time. As the interpolating operators couple to continuum states such as N π, it is possible that, by using Lüscher's method (and ideally, multiple spatial volumes), the phase-shifts for meson-baryon scattering can be extracted in channels with contributions from disconnected diagrams. In addition to excited single baryon states, and the continuum states that carry zero units of momentum in the volume, there are also continuum states where each hadron carries one or more units of momentum in the volume, while having vanishing total momentum. The lowest energy state containing hadrons A and B with back-to-back momenta ±p = ± 2π L n (where n is an integer triplet) occurs at In attempting to unravel the spectrum of states contributing to the correlation functions, we must also consider such continuum states. The lowest-lying negative parity state that is expected to couple to the interpolating operator for the single nucleon is the s-wave N π state (more precisely, we refer to the A + 1 representation of the hyper-cubic group), which has a threshold, neglecting interactions, of M π +M N = 0.27618±0.00034±0.00011. Fitting the EM shown in fig. 21 between time-slices t = 93 to t = 119 yields, significantly above threshold. Therefore, we conclude that this state is an s-wave πN scattering state with isospin I = 1 2 , as it has energy considerably below that of the first momentum excitation in the volume at E |n|=1 N π = 0.33889 ± 0.00042 ± 0.00014, and the ππN state. Given that there are no other channels that are energetically allowed for this state to mix with, the s-wave πN phase-shift 8 in this channel can be extracted at an energy of where the uncertainties are dominated by the uncertainty in E N ( 1 2 − ) . It is important to note that this channel receives contributions from disconnected diagrams, and in the calculation we are doing, these contributions are completely accounted for in the gauge-configurations. Using the standard Lüscher procedure, a phase-shift of δ πN = −26 ± 7 ± 6 degrees is found at this energy. The Prony histogram in fig. 21 shows significant structure in this channel, and one could argue that there is a single level at E ∼ 0.45, but this would require further exploration.
For the negative-parity state that couples to the interpolating operators for the Λ, the situation is not so clean. The thresholds for the lowest-lying continuum states, Σπ and N K, are located at M Σ + M π = 0.29747 ± 0.00030 ± 0.00019 and M N + M K = 0.30384 ± 0.00033 ± 0.00011, respectively, in the absence of interactions. The corresponding lowestlying states with one unit of momentum occur at E |n|=1 Σπ = 0.35857 ± 0.00037 ± 0.00023, and E |n|=1 N K = 0.35763 ± 0.00039 ± 0.00012. Fitting the EM shown in fig. 22 between time-slices t = 88 to t = 117 yields, This is, within uncertainties, at the threshold for Σπ or N K. The eigenstates will be a combination of these two systems and it is likely that we have not resolved the two nearbystates in the EM, and the result in Eq. (36) is actually an average of two closely-spaced energies.   For the lowest-lying negative-parity state(s) produced by the interpolating operator for the Σ, the situation is even more complicated. The threshold of the non-interacting swave Σπ state is at M Σ + M π = 0.29747 ± 0.00030 ± 0.00019, for the Λπ state is M Λ + M π = 0.29191 ± 0.00030 ± 0.00007, and for the N K state is M N + M K = 0.30384 ± 0.00033 ± 0.00011. Therefore, in this large volume, we expect to observe three eigenstates that are nearly degenerate. The lowest-lying states with one unit of momentum occur at E |n|=1 Σπ = 0.35857 ± 0.00037 ± 0.00023, E |n|=1 Λπ = 0.35341 ± 0.00030 ± 0.00007, and E |n|=1 N K = 0.35763 ± 0.00039 ± 0.00012 and are well separated from the n = 0 states. Fitting the EM shown in fig. 23 between time-slices t = 95 to t = 118 yields, in the region where one expects to find three closely-spaced states, corresponding to the eigenstates dominated by Σπ, Λπ, and N K. Given how closely spaced these states are expected to be, the extraction in Eq. (37) is likely a complicated average of three energies. The situation is no better for the lowest-lying negative-parity states that are expected to couple to the interpolating operator for the Ξ. The lowest-lying s-wave continuum states are πΞ, ΛK and ΣK. The threshold for these states, in the absence of interactions, are E = 0.37731 ± 0.00034 ± 0.00021, respectively, Therefore, we expect to observe two sets of three nearly degenerate eigenstates. Fitting the EM shown in fig. 24 between time-slices t = 91 to t = 118 yields, in the region where one expects to find three closely-spaced states, corresponding to the eigenstates dominated by n = 0 Ξπ, ΛK, and ΣK. Given how closely spaced these states are expected to be, the extraction in Eq. (38) is likely an average of three unresolved energies.
There are hints of a couple of other peaks in the Prony-histogram, but nothing conclusive. The results of the best extractions of the ground-state baryon masses using multiexponential fitting and the matrix-Prony method, which give consistent results for each species of baryon, are collected in Table VI. The masses of the lowest-lying J π = 1 2 − states with unit baryon number extracted by fitting three exponentials and by the matrix-Prony method. The first uncertainty is statistical while the second is the fitting systematic.

VIII. SIGNAL-TO-NOISE RATIOS
Many observables of importance to particle physics that are currently being calculated with LQCD, such as the pion decay constant and the Gasser-Leutwyler coefficients, require the calculation of mesonic correlation functions. Statistical fluctuations on each time-slice of these correlation function are well-behaved. In contrast, as argued by Lepage [10], correlation functions involving one or more baryons exhibit exponentially growing statistical noise. In the case of a single positive parity nucleon, the correlation function has the form where N is an interpolating field that has non-vanishing overlap with the nucleon and the angle brackets indicate statistical averaging over measurements on an ensemble of configurations. The variance of this correlation function is and therefore, as Lepage [10] argued, the noise-to-signal ratio behaves as More generally, for a system of A nucleons, the noise-to-signal ratio behaves as Therefore, in addition to the signal itself falling as G ∼ e −AM N t , the noise-to-signal associated with the correlation function grows exponentially, as in Eq. (42). These arguments are constructed for a system with an infinite time-direction and are modified in an important way for systems with a finite time-direction with given BCs. The calculations that are presented in this work have employed anti-periodic BC's in the time-direction. With such BCs the positive parity nucleon correlation function in Eq. (39) becomes where E N π is the energy of the lowest-lying negative-parity state in the volume, which, for this ensemble of configurations, is a continuum nucleon and pion at rest. The arrow denotes the behavior of the correlation function far from source (in both time-directions). Further, the correlation function dictating the behavior of the variance of the nucleon correlation function is modified similarly, with Eq. (40) becoming The first term in Eq. (44) arises from 3π's propagating forward and 3π's propagating backwards, the second term arises from 2π's propagating forward along with one π propagating backward and vice versa, the third (time-independent) term arises from a nucleon propagating forward and a nucleon propagating backward, and the ellipses denotes terms involving larger masses. As the negative-parity state is more massive than the nucleon, the nucleon is the dominate component in the correlation function, Eq. (43), for a number of time-slices beyond the mid-point of the configuration. From this argument, one expects to see the signal-to-noise ratio degrade even more rapidly than the expectation shown in Eq. (41) in the time-slices near the mid-point of the configuration where the correlation function is still dominated by the nucleon. One expects to find regions of the correlation function, depending on the structure of the source, which have the noise-to-signal scaling as MπT e (mp+ 3 2 Mπ)t , and e mp(t−T ) , or combinations thereof.
The high statistics calculations we are presenting here enable a detailed study of the behavior of the signal-to-noise ratio associated with the correlation functions formed with quark propagators generated with anti-periodic BCs.
It is useful to form the effective noise-to-signal plot, in analogy with the EMs. 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. 25, the full EM of the smeared-point nucleon correlation function is shown (with t J = 3), and in fig. 26, the full EM of the smeared-smeared nucleon correlation function is shown (with t J = 5). Also shown are the energy-scales associated with the growth of the noise-to-signal ratio from Eq. (46), with uncertainties generated using the Jackknife procedure. Considering the smeared-point correlation function in fig. 25, after time-slice t = 35 or so, the correlation function is dominated by the ground-state nucleon which persists until time-slice t ∼ 70. Beyond this time-slice the backward propagating negativeparity N π-state becomes dominant. Between time-slices t ∼ 40 and t ∼ 50, the noiseto-signal ratio is determined by the expectation of m p − 3 2 M π . However, after t ∼ 50 the signal-to-noise ratio degrades exponentially faster than this, and by t ∼ 65 the relevant energy-scale is ∼ m p + 1 2 M π and increasing with t. Similar behavior is clear in the smearedsmeared correlation function, for which the nucleon ground state dominates from an earlier time-slice.
It is clear from this analysis of the noise-to-signal ratio, that the length of the timedirection of these configurations and resulting thermal states are limiting the precision of the ground-state nucleon mass determination. This will be even more true for the multiple baryon correlation functions for which the signal-to-noise degrades exponentially faster than in the single nucleon correlation functions. Increasing the length of the time direction will lead to exponential improvement of the correlation function at large times where the nucleon component dominates the correlation function. It is interesting to note that the coefficients of the backward propagating contributions to the noise-to-signal ratio are suppressed by powers of e − 1 2 MπT . On the current configurations with T = 128, in order to reduce the contribution to the noise from the m p − 1 2 M π component by an order of magnitude, the time extent would need to be increased to T ∼ 192. This will reduce the m p + 1 2 M π component by a factor of ∼ 84 and the m p + 3 2 M π component by ∼ 770. Such an increase in the temporal extent would significantly decrease the statistical uncertainties with which ground-state signals are extracted.
The noise-to-signal analysis of the Ξ correlation functions is somewhat more complex, because there are a number of low-lying states which can contribute to the variance. For the ΞΞ noise correlation function, the lightest intermediate states that can couple to quarkcontent ssussu are KKη, ηηπ and ηηη. The full EMs and the E S plots for the smearedpoint and smeared-smeared Ξ correlation functions are shown in fig. 27 and fig. 28. In both the smeared-point and smeared-smeared Ξ correlation functions, the noise-to-signal ratio is growing exponentially slower than naive expectations, until about time-slice t ∼ 55. As the Ξ ground-state dominates the smeared-smeared correlation function beyond t ∼ 40, this allows for a extraction of the mass with higher precision than expected. This suggests that the noisesource does not couple to the low-lying mesonic states as strongly as expected, and that more massive mesonic states are dominating the noise over many time-slices. However, eventually, for t 55, the growth of noise overshoots the original Lepage expectation (indicated by the lowest horizontal line in figs. 27 and 28. Whilst we are primarily interested in noise in the baryonic sector, it is interesting to note that the mesonic correlation functions also suffer from similar issues. According to the above arguments, the pion correlation function on lattices at zero temperature (infinite temporal extent) will have noise that is independent of time (up to fluctuations) while the kaon will have noise that grows exponentially with the small energy difference m K − 1 2 m π − 1 2 m η . However at finite temperature, the noise correlation functions of both systems receive additional contributions that grow faster than the above expectations. This is shown in fig. 29.

IX. SCALING WITH COMPUTATIONAL RESOURCES
An important component of our current work is to address the future requirements for LQCD calculations in nuclear physics, a field characterized by small energy scales in heavy systems, for example, the 2 MeV binding energy of the ∼ 2 GeV deuteron. In fig. 30, we show the extracted mass of the π + , K + , N and Ξ as a function of the number of configurations in the ensemble for both the exponential and matrix-Prony analysis methods. The full set of measurements performed on each configuration are included, and the fitting intervals are chosen to optimize the extraction for each ensemble size. In each case, the uncertainty in the mass is reduced, as expected, with increasing ensemble size, and the mass extracted from the smaller ensembles tends to be less than that from the larger ensembles. Fig. 31 shows the fractional uncertainty in the mass of the π + , K + , N and Ξ, associated with the results in fig. 30, as m Ξ , and m Ξ + m η − 1 2 M π (from lowest energy to highest energy).
The dependence of our results for hadron masses on the number of sources used in the calculations is explored in fig. 32 where we show the fractional uncertainty in the mass of the π + , K + , N and Ξ as a function of the number of sources used on each configuration. In this figure we use an ensemble of 1012 configurations, those on which there are at least 100 measurements. A simple fit of the form δM/M = AN b src returns exponents b = −0.03(2), -0.65 (19), -0.41(3), -0.40(6) for the π + , K + , N and Ξ, respectively.
The results of this analysis can be simply summarized for baryons (averaging over the nucleon and Ξ) as For mesons, a similar scaling is seen, with a somewhat worse scaling with the number of measurements per configuration in the case of the pion, consistent with the saturation seen in fig. 3. This functional form enables us to quantify the relative benefit of increasing the number of sources per configuration compared to increasing the total number of configurations. The costs involved in this are as follows: • Gauge configuration generation: The total cost of generating the ensemble of 1194 gauge configurations was 2M JLab-6n cluster node-hours and the production took a significant amount of wall-clock time. Configuration generation costs scale linearly with the number of configurations once a Monte-Carlo trajectory has thermalised (in this case the overhead of thermalisation was approximately 10%). In order to generate significantly larger ensembles (containing 10 4 or 10 5 gauge fields) in a reasonable wallclock time, it will be necessary to run multiple trajectories in parallel. Given wall-clock time and memory constraints, an individual trajectory will produce O(1000) gaugefield configurations that are useful for measurements. Consequently the thermalization overhead will conservatively remain at about 10%. Each configuration requires ∼ 2 × 10 3 JLab-6n node-hours to produce.
• Measurement calculations: The total cost of computing all of the measurements performed in this work was 7 × 10 6 Jlab-6n node-hours. The cost to generate the 245 light-quark and strange-quark propagators per configuration on the 1194 configurations in this ensemble was ∼ 3M Jlab-6n node-hours, while the cost to generate the baryon and meson blocks (used at intermediate stages of the calculations) was ∼ 3.5M Jlab-6n node-hours. Contracting the blocks to accomplish the desired measurements (one, two, ... baryons, one, two,... mesons and so forth) cost ∼ 0.5M Jlab-6n nodehours. If propagators on a given configuration are computed in sets of 100, the initial overhead of constructing deflation vectors in the EigCG algorithm becomes negligible (at the 1% level) and can be eliminated for further sets of calculations by storing the eigenvector information. On typical machines, each set of propagators and associated hadron blocks (technically, not an efficient way to calculated the single hadron spectrum, but critical for two and more hadron calculations) requires 22 JLab-6n node-hours to produce.
• Anisotropy: The anisotropy of the lattices used in our calculations proved useful in reducing systematic errors in our fits (see Table III), providing approximately a 1/ ξ ef f reduction. However, the cost of producing gauge-field configurations and propagators scales as approximately ξ 2 for the same physical extent (one power arises from the additional time-slices and one power arises from the worsening condition number of the Dirac operator). Comparing these exponents, we would conclude that using anisotropic configurations is not ideal. However for more complicated multi-hadron systems where useful fit ranges are much reduced in physical units, the anisotropy will likely prove to be very useful. This remains to be investigated further in subsequent studies.
Using this information and the scalings in Eq. (47), we can address the question of how much computation is required to achieve a particular level of statistical precision. With the current data this is only possible in the single hadron sector; ongoing analyses will address the B > 1 in the near future. To halve the uncertainty in the determinations of the ground- state baryon masses (calculating the nucleon mass at the ∼ 1 MeV-level), an increase in the number of configurations by a factor of four, or of the number of measurements per configuration by a factor of 5.6 is required. Achieving this precision by performing more measurements on the existing set of configurations (∼ 1100 additional measurements on 1200 configurations) would cost 30 M JLab-6n node-hours. Achieving the same precision by generating an additional configurations and performing the same number of measurements on them (3600 configurations with 245 measurements) would require 27 M JLab-6n node-hours. Both approaches have similar cost at this level of precision, but for further improvements, the generation of additional configurations will be more efficient. Additionally, the second approach will further improve the uncertainty for observables such as the pion mass that have saturated in terms of the number of sources per configuration. This approach would clearly be of more benefit to the broader community.

X. CONCLUSIONS
The energy-scales that arise in nuclear physics are typically in the MeV range, and in order for LQCD to have significant impact in this field, baryon masses (and energy-eigenstates in the volumes relevant to scattering processes) must be calculable with uncertainties that are a fraction of an MeV (including isospin-breaking and electromagnetic interactions, quark mass, lattice volume and lattice spacing extrapolations). Current computational resources do not permit such calculations. In this work we have performed the first high-statistics study of baryon correlation functions to better understand a number of issues that will impact the precision with which quantities of importance to nuclear physics can be determined with LQCD. In the future, we will extend our analysis to look at observables in the B > 1 baryon sector.
At the single lattice spacing, lattice volume and unphysical light quark mass used in this work, we find the following set of ground state masses which we present in physical units. Since the lattice spacing is known with less precision than the lattice masses presented here, we make the systematic uncertainty arising from the lattice spacing explicit (third uncertainty). Given that there are significant ambiguities in scale setting, the most precise result will be for dimensionless quantities.
With high precision measurements of baryon correlation functions obtained from a single type of source for the light-quark and strange-quarks propagators, we have shown that the number of methods that can be used to extract the arguments of the contributing exponentials increases. This is due to the fact that some methods become stable when the uncertainties become small, such as the method of Prony and also the direct fitting of multiple exponentials. Histograms constructed from the roots found in the Prony method are found to be useful in identifying mass regions where states may exist, but we have not yet arrived at a well-defined (rigorous) method with which to use these histograms directly. It is likely that a more refined statistical analysis of these correlation functions using the most modern statistical tools will further increase the physics that can be extracted.
The exponential signal-to-noise degradation that plagues baryon correlation functions is currently a serious limitation for the calculation of nuclear physics observables, and is one significant difference between particle and nuclear physics LQCD calculations. The high statistics calculations we have performed have allowed us to systematically explore this issue. We find that the issue is more serious than one would naively expect, due to (what in hindsight is now obvious) the use of anti-periodic BC's in the time-direction on the quark-propagators. The variance of the correlation function is symmetric about the mid-point of the time-direction of the configuration. Therefore, the optimal region in which to determine the baryon masses (and also their interactions) is in the first half of the configuration, far from the midpoint. This significantly reduces the number of useful time-slices. Given that most the time required for these calculations is in the measurements, and not in the configuration generation, a cure for this problem is to generate ensembles of configurations that are longer in the time-direction than those currently being used (as opposed to working with different BC's on the quark propagators that are less theoretically "clean"). 9 The multi-exponential fitting and Prony methods enable the ground-state to be probed closer to the source where the statistical uncertainties are exponentially smaller (also one of the important aspects of the variational method), somewhat reducing the impact of the exponentially degrading signal-to-noise near the mid-point of the configuration. Given that the signal-to-noise degradation is exponentially more severe in systems containing two or more baryons, all currently available tools will be required to make optimal use of the computational resources. For excited states in a given channel, variational methods seem to be superior to the standard approach used here.
An important result of this work has been to quantify the statistical scaling of simple observables in the sub-percent regime of uncertainty. Scaling with the number of configurations was found to adhere to the expected 1/ N cfg behaviour. We have also investigated the issue of saturation, asking many measurements can be performed on a single gauge-field configuration before it becomes more cost effective to generate another statistically independent gauge-field configuration? To address this we have looked for deviations from 1/ √ N src behavior in the uncertainties in the correlation functions and extracted masses as a function of the number of measurements performed on each configuration in the ensemble. The measurements of the mesons start to saturate after a relatively small number of measurements, in this case of order ∼ 10, while the baryon correlation functions show no signs of saturation up to ∼ 200.
A natural question to ask is if the same extractions could have been performed with fewer computational resources by using the variational method with a number of different sources for the baryons. We estimate that comparable resources would have been required to achieve comparable uncertainties in the states we have examined. However, we have not been able to extract excited states of the nucleon with much precision because of closely spaced states with the same quantum numbers. This is likely to be something that the variational method would better control. Given that the present work was exploratory in nature, this is not a concern at present, but it is clear that high-statistics calculations of correlation functions arising from multiple interpolating operators will be required in order to explore the structure and interactions of nuclei. We are working on implementing this, but it will require significant computational resources to perform, even at pion mass ∼ 390 MeV and for a relatively small lattice volume and relatively coarse lattice spacing.
It is clear that sub-MeV uncertainties an hadron energies will become routine with the anticipated increase in computational resources available to lattice QCD, and that the small energy scales that characterize nuclear physics are within reach. However, this program will require large ensembles of gauge-field configurations that have large extent in the timedirection, and will require a large fraction of the computational resources devoted to measurements.

XI. ACKNOWLEDGMENTS
We thank R. Edwards and B. Joo for help with the QDP++/Chroma programming environment [13] with which the calculations discussed here were performed. KO would like to thank A. Stathopoulos useful discussion on numerical linear algebra issues and for his contribution in the development of the EigCG algorithm. 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