Evolutionary Fitting Methods for the Extraction of Mass Spectra in Lattice Field Theory

One of the most important and common uses of curve fitting in lattice gauge theory is the extraction of hadronic masses and matrix elements from measured correlation functions.

The problem of extracting the mass spectrum from the correlators measured in Lattice QCD simulations is well known [1,2]. The simulation produces data G i = G(t i ) for the expectation values of the correlator G(t) at a finite number of discrete (Euclidean) time steps t i , 1 ≤ i ≤ N. On the other hand, from theory the exact form of the propagator is known to be given by an infinite series2

where we will assume that the energy levels are ordered, E n ≤ E n+1 . The problem is then to determine an infinite number of amplitudes Z n > 0 and energies E n > 0 from only a finite number of data points G i , an obviously ill-posed problem.

To make the problem well-posed, we have to add some further physical constraints. The piece of theoretical information that is normally used is that the sequence of the Z n is bounded from above, and therefore the correlator will be dominated by the lowest few terms at all but the smallest values of t due to the exponential suppression by E n . We can therefore truncate the series in equation (1) after a finite number of terms, provided we only attempt to fit at t > t min large enough for the truncation to include all terms that make a significant contribution.

In doing so, we are faced with a choice: Either take a fixed number n max of terms of the sum in (1) and adjust t min so as to obtain a good fit, or fix t min and vary n max so as to extract the largest possible amount of significant information. The problem with the first strategy is that we are throwing away valuable information contained in the data points at t < t min , leading to large statistical errors if t min is chosen too big. These need to be balanced against large systematic errors arising if t min is chosen too small for the given n max [1].

In this paper we will therefore adopt the second approach and attempt to fit all data points (excluding only the single point at t = 0 for practical reasons)

with a variable number n max of exponentials.

Naively, one might want to try to simply run a series of fits using a state-of-the art fitting method such as Levenberg-Marquardt [3,4] at a number of different n max , and with a variety of initial parameter values, and choose the fit that produces the lowest χ 2 per degree of freedom. This method can be made to work in the case of one single correlator if the problem of finding an appropriate starting point in a potentially multi-modal landscape is somehow solved.

However, for many questions in lattice QCD it is necessary to fit multiple correlators, which may or may not share some of the energy levels E n . In this case, the fast combinatorial growth with n max of the number of possibilities of assigning shared or separate fit parameters E n to the fitting functions for the different correlators renders the application of this naive method to those problems largely impossible. The problem of choosing appropriate initial parameter values further exacerbates this approach.

The current state of the art using the variable-n max approach is that taken in [1] where it is used in the context of a Bayesian method [5] of constrained curve fitting. The latter works by adding prior information about the fit parameters and using it to constrain the fit to a smaller subset of likely parameter values out of the space of all possible values. Only those parameters whose fitted values are largely independent of the priors used to constrain them are considered to be determined by the data, while the others are disregarded as having been imposed by the priors.

While it is thus possible to determine which fitted quantities are independent of, or only weakly dependent on, the priors and thus determined from the data, the idea of using some external knowledge as an input could be seen as incompatible with the notion of a first-principles determination of the quantities of interest from QCD itself, without any recourse whatsoever to empirical models. Moreover, in some cases appropriate priors may not be available. Under those circumstances, it becomes desirable to be able to extract some estimate of the parameters to be fitted from the data themselves, and to use this estimate as a prior in the context of a Bayesian constrained fit. A number of methods to do this, such as the sequential empirical Bayes method [6], have been used in the existing literature.

A completely different state-of-the-art approach that does not rely on prior information while allowing for extraction of a spectrum from multiple sets of data for improved statistics is the variational method [7,8,9]. Here one sets out by fixing a channel corresponding to a specific set of quantum numbers, and finds a set of appropriate linearly independent operators O i for th

…(Full text truncated)…