# A novel analysis method for excited states in lattice QCD - the nucleon case

###### Abstract

We employ a novel method to analyze Euclidean correlation functions entering the calculation of hadron energies in lattice QCD. The method is based on the sampling of all possible solutions allowed by the spectral decomposition of the hadron correlators. We demonstrate the applicability of the method by studying the nucleon excited states in the positive and negative parity channels over a pion mass range of about 400 MeV to 150 MeV. The results are compared to the standard variational approach routinely used to study excited states within lattice QCD. The main advantage of our new approach is its ability to unambiguously determine all excited states for which the Euclidean time correlation function is sensitive on.

## I Introduction

The study of excited states within the framework of Quantum Chromodynamics on the lattice (LQCD) is difficult since it is based on the evaluation of Euclidean correlation functions for which the excited states are exponentially suppressed as compared to the ground state. The standard approach to study excited states is based on the variational principle, where one considers a number of interpolating fields as a variational basis and defines a generalized eigenvalue problem (GEVP) using the correlation matrix computed within the chosen variational basis. The GEVP has been widely applied in the study of hadron spectroscopy by a number of lattice groups with recent results given in Refs. Basak et al. (2007); Gattringer et al. (2008, 2009); Mahbub et al. (2010); Bulava et al. (2010); Alexandrou et al. (2013a).

In this paper, we examine in depth the application of a new method based on statistical concepts for extracting the excited states from Euclidean correlators to study the nucleon spectrum. The so called Athens Model Independent Analysis Scheme (AMIAS) Papanicolas and Stiliaris (2012) was originally developed to extract scattering amplitudes from experimental measurements in the to transition Papanicolas and Stiliaris (2007) without utilizing a specific model thus avoiding model dependencies of the outcome. A first analysis of the nucleon two-point function in LQCD was carried out with promising results in Ref. Alexandrou et al. (2008). In this work, we extend the method to analyze correlation matrices and compare the results with those obtained in our recent study using the standard variational method Alexandrou et al. (2013a). An advantage of AMIAS is that it can be applied to the correlation function at small separation times by allowing any number of states to contribute rather than to the large time-limit behavior typically done in the variational method. In fact, the merit of the method is that it determines the actual number of states on which the correlation matrix is sensitive on. Thus AMIAS does not rely on plateau identification of effective masses, which are usually noisy and thus difficult to determine, but instead it utilizes all the information encoded in the correlation function with the advantage of exploiting the small time separations where the statistical errors are small.

The paper is organized as follows: In section II we give the general description of the method, including a discussion of an importance sampling version of the algorithm. The method is compared to standard variational approach paying particular attention to the determination of the number and correlations of the fitting parameter, which are advantages of our method. Section III gives the lattice QCD formulation pertinent to the problem at hand. Finally, in section IV, we discuss the results extracted using AMIAS and in section V we summarize our findings and give our conclusions.

## Ii The AMIAS approach

In this section we introduce AMIAS and demonstrate its applicability in extracting the excited states from Euclidean correlation functions.

The starting point of AMIAS is a set of measurements with the associated statistical errors Papanicolas and Stiliaris (2012). Since the quantity of interest in the work is a Euclidean two-point function computed at discrete times we explain the method based on the measurement of these correlators. Let us denote by the measurement at time . The next ingredient for AMIAS is a theoretical parameterization of the measurements or an ’underlining theory’. For lattice QCD correlators this is well-known from the spectral decomposition of the hadron propagator.

(1) |

The exact correlator is approximated as the average over measurements of an appropriately defined quantity to be described in Section III:

(2) |

for a total of measurements. The parameters and are to be determined by AMIAS, where we order the exponentials by the value of . For large values of the time the exponential with the smallest exponent dominates.

More details on the lattice formulation will be given in Section III.

The Central Limit Theorem states that the probability for the above average at time to have a value equal to is

(3) |

where is the standard deviation at time and is, for our purpose, an irrelevant normalization constant. The above result can be found in many textbooks and the assumptions made are that the measurements are uncorrelated and that is large. The value of in Eq. (3) will, in general, be underestimated if correlations are present. For LQCD the correct value of is usually obtained using either a jackknife or bootstrap procedure. An example of the distribution given in Eq. (3) for a given value of is shown in Fig. 1 showing indeed that the distribution is to a good accuracy a Gaussian. The jackknife error is found to be in agreement with that predicted by the Central Limit Theorem, which means that the measurements used for the average are uncorrelated.

The probability of obtaining a set of measurements, each for a different value of , is equal to the product of the probabilities for each data point:

(4) |

where we have set . If we use the expansion in terms of exponentials given in Eq. (1), then the above result can be written as

(5) | |||||

For practical purposes the infinite summation over the energy spectrum involving the parameters and must be truncated keeping a maximum number of energy levels:

(6) |

The appearing in the above equation is the usual used in the least-squares fitting method. The Probability Density Function (PDF) of Eq. (5) is the basis of the minimization in the case where the values of the parameters are unknown: Maximizing the probability is equivalent to minimizing the exponent.

The principal idea behind AMIAS Papanicolas and Stiliaris (2012) is that given a particular set of sampled averages, , any arbitrary value assigned to the parameters and constitutes a solution of Eq. 1 having the probability given by Eq. (5). Following standard statistical concepts Papanicolas and Stiliaris (2012), the probability that the parameter assumes a specific value in the range is equal to

(7) |

The above formulation has been applied to the analysis of experimental nuclear extrapolation data Papanicolas and Stiliaris (2007) and to a lattice QCD calculation of the nucleon excited states Alexandrou et al. (2008), using a uniform sampling for values of in Eq. (7), yielding consistent results with the standard analysis.

For reasons that will become apparent later on, instead of using a single correlation function we can apply AMIAS to a correlation matrix with the spectral decomposition

As in the case of a single correlator we keep a maximum number of energy levels in the sum.

The method can be applied to a general matrix and in the analysis where we apply AMIAS we will be interested in correlation matrices that are hermitian. This constrains the number of amplitudes , simplifying the fitting problem. Although the matrix elements can be correlated with each other Eq. (6) is valid for each matrix element independently, hence the PDF for the entire correlation matrix can be expressed as the product of the PDF for each matrix-element, leading to a that is the sum of the individual ’s:

(8) |

### ii.1 Importance sampling

While uniform sampling works well for a single correlator, for a correlation matrix determining the fit parameters and is impractical due to the large number of parameters. One can use the Metropolis algorithm Metropolis et al. (1953) to sample more efficiently the distribution , through a random walk. However, the Metropolis algorithm being a local random walk has the disadvantage of getting trapped at local minima. In order to avoid this problem and ensure access to the full multidimensional space, we employ parallel tempering Swendsen and Wang (1986).

In our parallel tempering scheme we run a number of additional
random walks, randomly initialized at different ‘temperatures’. Thus we have a chain of
random walks sampling the sequence of PDFs:

,
where the s are the parallel tempering temperatures. Sampling the
distribution histogram of each parameter is done through
Eq. 7 in the usual way but only the original PDF is finally used
obtained by setting . The high temperature PDFs are generally able to
sample larger phase space, whereas the exact PDF with
whilst having precise sampling in a local region of phase space, may
become trapped in local minima. Information from the high temperature walks is passed down to the low
temperature ones through exchanges among ensembles. Swaps are
normally attempted between systems with adjacent temperatures, and
, and are accepted with probability

(9) |

Parallel tempering is an exact method, in that it satisfies the detailed balance condition.

### ii.2 Number of parameters

One of the significant advantages of AMIAS is that it determines unambiguously the parameters to which the data are sensitive on i.e. it determines in the truncation of the infinite sum in Eq. (6). The results obtained are then invariant under changes of . The strategy is to increase until there is no sensitivity to the additional exponentials and thus no observable change in the sampled spectrum. In Fig. 2 we show an example of such an analysis for the nucleon correlator , which will be defined below. As we increase , additional exponential terms are identified resulting in a well defined Gaussian-like distribution for each of the parameters. We can get values for the parameters, , and , by fitting the distribution of each mass to a Gaussian. For exponential terms beyond the first three (e.g. those having exponents and ) we get a uniform distribution indicating that there is no contribution to the minimization of from these terms. We can thus safely conclude that is a safe choice since it includes, in addition to the well determined exponentials, an additional insensitive exponential term. Table 1 gives the numerical values of the energy eigenstates corresponding to the distributions of Fig. 2, where we also compare with a standard least squares minimization algorithm. The advantage of AMIAS is in the consistency of the results even when exponential terms are included on which the data are insensitive. With the usual least squares minimization one can identify three states with an error of over 50% for the . No indication as to the presence of higher states can be extracted. AMIAS on the other hand yields with a 10% error. As we increase the number of exponentials from 3 to 5 the results for the parameters on which the correlator has sensitivity remain unchanged demonstrating the robustness of the method.

AMIAS | Standard Least Squares | |||||
---|---|---|---|---|---|---|

2 | – | – | ||||

3 | ||||||

4 | – | – | – | |||

5 | – | – | – |

### ii.3 Fit range

Another parameter that enters in the analysis of correlators is the initial time used in the fits to the exponential form. The larger the initial time the smaller is the contribution of exponentials beyond the first one. Thus, we can suppress the contribution from the higher exponents by varying the starting value . In the standard analysis the time is varied until the so-called effective energy defined by

(10) |

becomes a constant as a function of (plateau region) yielding the lowest exponent which will be identified as the ground state energy in our lattice QCD study. Similarly, in the AMIAS analysis we need to verify that the determination of the exponents of the first exponential terms will be unaffected as we vary . In Figs. 4 and 4 we show such an analysis for both and is performed. By increasing the three exponential with exponent is eliminated but the previous two exponents remain unaffected. Elimination of the third exponent is also apparent when examining the distribution of the amplitudes where the distribution of becomes uniform as is increased. Therefore, by increasing one removes unnecessary correlations from terms that we cannot clearly obtained within the available statistics and in addition, one verifies the validity of the sampled dominant exponentials.

### ii.4 Correlations

A central issue that is properly treated in AMIAS, is the handling of correlations, since all possible correlations are accounted for. The sampling method of Eq. (7) allows all fit parameters to randomly vary and to yield solutions with all allowed values, including the insensitive exponential terms. The visualization of the dominant correlations can be accomplished by a two-dimensional contour plot in which the values of a selected pair of parameters is varied around the region of maximum probability, keeping all other parameters fixed. In Fig. 5 we present such a correlation analysis, where the plane defined by the values of the parameters is color coded according to the -value. The top-left and bottom-left parts are examples where the parameters are correlated, while the right part is an example of uncorrelated parameters. In particular the bottom-right part is an example where sampling is insensitive to one of the parameters, namely , indicating in this case the absence of a fourth state (compared with Fig. 2).

## Iii Lattice techniques

In this section we give a more detailed description of the lattice QCD results that will be analyzed using AMIAS.

### iii.1 Simulation details

We use gauge configurations produced using the twisted mass fermion action with two degenerate light quarks (). We also use twisted mass fermion gauge configurations adding to the light quarks a strange and charm quark with masses fixed to their physical values (). More details on the gauge configurations produced by the European Twisted Mass Collaboration (ETMC) can be found in Ref. Baron et al. (2010) and for the in Ref. Boucaud et al. (2008). The twisted mass fermion ensembles used in this work are summarized in Table 2. The pion mass values span a range from about 210 MeV to 450 MeV. Apart from the twisted mass fermion ensembles we also analyze an ensemble of clover fermion gauge configurations produced by the QCDSF collaboration with near-physical pion mass of MeV Bali et al. (2013). The nucleon excited states have been analyzed on these gauge configurations using the conventional approach in Ref. Alexandrou et al. (2013a) and thus serve as a good set for making detailed comparison with the results obtained with AMIAS.

ensembles | |||
---|---|---|---|

Twisted mass fermions, fm. | |||

, | (GeV) | 0.2696(9) | 0.3082(6) |

fm | No. of confs | 659 | 232 |

clover fermions, fm. | |||

, | (GeV) | 0.160 | |

fm | No. of confs | 250 | |

ensembles, fm | |||

, | (GeV) | 0.375162 | |

(B55.32) | No. of confs | 5000 |

### iii.2 Correlator functions

In order to study the energy spectrum within the framework of lattice QCD one evaluates the Euclidean two-point correlation function

(11) | |||||

where is a creation operator acting on the QCD vacuum (interpolating field) having the same quantum numbers as the states of interest and is the three-momentum. The two-point correlation function can be expressed as a sum of the energy eigenstates of QCD that exponentially decay as a function of time with an exponent that depends on the energy of the state. Thus by fitting to such a form one can in principle extract the energies of the system. The problem, however, is that the higher excited states decrease exponentially as compared to the ground state and extracting them is difficult since the signal to noise decreases exponentially for all hadrons except the lowest mass pseudoscalar. Fitting such correlation functions to two-exponentials can be one approach to extract the first excited state when one has lattice data with small enough errors. On the other hand, extracting the ground state is much easier since asymptotically it is the only state that survives the large time limit of .

Let us consider the zero momentum correlator obtained from Eq. (1) by setting . In a simulation we have measurements at each time , with the measurement having the form , where is the number of lattice time slides on which the correlator is evaluated. Note that from now on we drop the discrete index on . The measurements are over a representative ensemble of gauge configurations generated via a Monte-Carlo sampling of the probability density function of the Euclidean lattice QCD action. For each value of we can form the average , which was the example used in the previous section to illustrate AMIAS.

The results shown in Figs. 4, 4 and 5 are carried out using the B55.32 ensemble of twisted mass fermions of a pion mass MeV, for which a large number of gauge configurations are generated. We will use this ensemble to perform a detailed comparison of AMIAS with the standard analysis used in lattice QCD for the study of excited states.

The standard method to study excited states in lattice QCD is the variational approach. One expands the basis to interpolating fields of the quantum numbers of states of interest and construct a correlation matrix. In this work we consider a correlation matrix of the form

where we denote the different types of interpolating fields with two type of indices and . The superscripts on the correlation matrix correspond to different levels of Gaussian smearing, while the subscripts to different spin combinations of the nucleon interpolating fields given by:

(12) |

Gaussian smearing is applied to the quark fields using the hopping matrix :

(13) | |||||

We also apply APE-smearing to the gauge fields entering . The parameters for the Gaussian smearing and are optimized using the nucleon ground state Alexandrou et al. (2013b). The local nucleon interpolator, , is well known to have a good overlap with the ground state of the nucleon. The trace in Eq. (III.2) is taken over Dirac indices and the correlation matrices and yield the positive and negative parity states of the nucleon, respectively. The states are eigenstates of the Hamiltonian with and we assume that the temporal extent of the lattice is large enough to neglect boundary contributions.

The low-lying energy spectrum should be unaffected by what nucleon interpolating field is used and in particular how much smearing is applied as determined by the parameters and . This is because Gaussian smearing only effects the coupling of the interpolating field with the nucleon eigenstates i.e. only the are affected but not the . As long as the amplitudes are non-zero the energy of interest can be extracted. In Fig. 6 we verify the independency of the energies by comparing the results from a smeared-smeared correlator where both interpolating fields in Eq. (II) are smeared to results extracted from a local-smeared correlator , where only the interpolating field at the source is smeared. We show results for the ground and first excited states. As expected the distributions for and are statistically equivalent but the local-smeared correlator provides a better estimate for the excited state since it has less overlap with the ground state. In contrast, the amplitudes and depend on the level of smearing. This is an indication that reducing smearing will improve the results extracted for the excited states, since their contribution to the correlator will be larger.

### iii.3 Comparison of AMIAS with the variational method

The standard way to extract the ground state energy in lattice QCD is to probe the long time-limit of the correlation function and consider the effective energy , which becomes independent of when the ground state is the dominant contribution to the correlator. Fitting the effective energy (or mass) to a constant in the plateau region yields the ground state energy (or lowest mass). As has already been demonstrated in Ref. Alexandrou et al. (2008), the value obtained with AMIAS for the ground state is in agreement with the one extracted from the effective mass plateaus. In this section, we focus our attention to the extraction of the excited states.

In order to study excited states in lattice QCD, one usually applies the variational method and defines a generalized eigenvalue problem (GEVP) given by

(14) |

where . The corrections to decrease exponentially like where Lüscher and Wolff (1990) for fixed . In a recent application of the variational method for the analysis of the nucleon states Alexandrou et al. (2013a) we found that a variational basis constructed from different smearing levels of the standard interpolating field containing both a small and a large number of Gaussian smearings is an appropriate basis for extracting the first positive parity excited state, known as the Roper.

In order to compare AMIAS with the results extracted using GEVP, we consider the same variational basis used in that study, namely we set the values of the smearing parameters to and 10, 30, 50, 180 and 300. These different smearing levels are labeled by the superscript on and produce a source with a root mean square radius and in units of the lattice spacing , respectively. The resulting matrices are symmetrized. For this analysis, we use 200 twisted mass configurations produced at =3.9, a = 0.004 or m MeV on a 32 64 lattice. In Fig. 7 we show the results extracted from AMIAS for the various correlation matrices. The value for the first excited state obtained from fitting the diagonal elements is statistically equivalent for all values of the smearings. This is illustrated in Fig. 7, where we show the results extracted from and , which correspond to the smallest and largest smearings, respectively. When a correlation matrix is used there is a lowering in the value of the first excited state and in addition a second excited state appears close-by. This indicates that this state cannot be detected in the diagonal correlators, presumably having a very small overlap within the standard nucleon interpolating field. A correlation matrix using the first three smearing levels is compared with a matrix, . Using a variational basis that includes small and large smearings i.e and results in a lowering of the energy of the first excited state. Furthermore, using the full correlation matrix, , results in a reduction of the error but leaves the mean values unaffected.

We compare the results extracted with AMIAS to the effective energies extracted from GEVP in Fig. 8. The AMIAS values are shown by the bands, while the effective mass by the symbols. As can be seen, the effective energies as a function of the time separation become very noisy for the first and second excited states making the identification of the plateau region ambiguous. AMIAS nicely extracts these energies, which are in agreement with the levels that the variational approach seems to indicate. This illustrates the advantage of AMIAS in the determination of the excited states.

## Iv The low-lying nucleon spectrum

Having presented a detailed comparison with the results obtained by using GEVP, in this section we present the results for the nucleon excited states in the positive- and negative-parity channels obtained by applying AMIAS. For this analysis we use the correlation matrix . These correlation matrices were also utilized within the GEVP analysis of Ref. Alexandrou et al. (2013a) and thus readily provide a comparison with the results obtained within the AMIAS analysis.

In Figs. 9 and 11 we show the results for both positive and negative parity states, for the twisted mass ensembles and for the clover ensemble analyzed in this work. The results are compared with those obtained using the GEVP. For the case of the clover action and the heaviest of the twisted mass ensembles (GeV) we have refined our analysis by expanding the basis of the GEVP as compared to the analysis carried out in Ref. Alexandrou et al. (2013a). For these two ensembles the results obtained with AMIAS are consistent with the GEVP. For the twisted mass ensemble with the smaller pion (GeV) the first excited state obtained with AMIAS has a larger error than the one extracted from the plateau. This is because in the GEVP one needs to fit for small in order to have enough fitting points. Contrary in AMIAS we utilize the whole time dependence and thus we are more confident for the AMIAS determination.

For the negative parity channel there is an interesting feature connected to the scattering state. Unlike the behavior observed in the positive channel where AMIAS yields results that did not show any dependence on here this is no longer the case. The lowest level extracted when using is higher that the energy of the S-wave scattering state assuming that no interactions take place. AMIAS yields results in agreement with the GEVP only after increasing the value of , albeit with larger errors. This was further investigated in Fig. 10, where we show the effective mass for the negative parity ground state for the three lowest pion masses used. The values extracted using AMIAS for various starting values of the fit range, are shown by the shaded error bands as . The values of used are indicated by the timeslice where the bands start. The dashed line shows the energy of the scattering state being the sum of the pion and nucleon mass. As can be see, with increasing the error band width increases becoming consistent with the mass of the scattering state but at the same remains consistent with the band obtained with . This is consistent with our analysis in section II.3 where we have shown that increasing suppresses excited states but does not affect the mean value of the ground state. One possible reason that one needs to increase to obtain agreement with the energy maybe the small overlap this two-particle state may have with a single particle interpolating field, requiring the suppression of the excited states to get a signal and may explain why other groups have not been able to detect it. In Fig. 11 we give the values extracted for the ground and the first excited state in the negative channel. The values shown for the correspond to the magenta bands in Fig. 10 obtained using for clover fermions and for twisted mass fermions. In contrast, the excited state was obtained using in analogy with the positive parity channel.

## V Summary and Conclusions

A novel method for the analysis of excited states within lattice QCD is applied to study the spectrum of the nucleon. The method uses importance sampling to probe the correct minimum of the multidimensional space defined by the states of the correlation matrix. AMIAS determines the number of excited states, which can be extracted from the information that the lattice data encode. It takes advantage of the whole time dependence of the correlators not requiring identification of the large time asymptot that limits the accuracy of the determination of excited states. The method is applied to successfully extract the first excited states in the positive and parity nucleon channels, identifying the Roper in the positive parity channel. In the negative parity it can extract the energy of the first excited state using the entire time dependence of the two-point function the energy of the scattering state is only extracted after eliminating the initial few time slices suppressing the excited states. This is a feature that we plan to explore further in future studies.

## Acknowledgments

We thank G. Koutsou for providing the two point functions for the clover fermion ensemble. Numerical calculations have used the Cy-Tera facility of the Cyprus Institute under the project Cy-Tera (NEA YOOMH/TPATH/0308/31) funded by the Cyprus Research Promotion Foundation.

## References

- Basak et al. (2007) S. Basak, R. Edwards, G. Fleming, K. Juge, A. Lichtl, et al., Phys.Rev. D76, 074504 (2007), eprint 0709.0008.
- Gattringer et al. (2008) C. Gattringer, L. Y. Glozman, C. Lang, D. Mohler, and S. Prelovsek, Phys.Rev. D78, 034501 (2008), eprint 0802.2020.
- Gattringer et al. (2009) C. Gattringer, C. Hagen, C. Lang, M. Limmer, D. Mohler, et al., Phys.Rev. D79, 054501 (2009), eprint 0812.1681.
- Mahbub et al. (2010) M. Mahbub, A. O. Cais, W. Kamleh, D. B. Leinweber, and A. G. Williams, Phys.Rev. D82, 094504 (2010), eprint 1004.5455.
- Bulava et al. (2010) J. Bulava, R. Edwards, E. Engelson, B. Joo, H.-W. Lin, et al., Phys.Rev. D82, 014507 (2010), eprint 1004.5072.
- Alexandrou et al. (2013a) C. Alexandrou, T. Korzec, G. Koutsou, and T. Leontiou (2013a), eprint 1302.4410.
- Papanicolas and Stiliaris (2012) C. N. Papanicolas and E. Stiliaris (2012), eprint arXiv:1205.6505.
- Papanicolas and Stiliaris (2007) C. N. Papanicolas and E. Stiliaris, AIP Conf. Proc. 904, 257 (2007).
- Alexandrou et al. (2008) C. Alexandrou, C. N. Papanicolas, and E. Stiliaris, PoS LATTICE2008, 099 (2008), eprint 0810.3982.
- Metropolis et al. (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
- Swendsen and Wang (1986) R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986), URL http://link.aps.org/doi/10.1103/PhysRevLett.57.2607.
- Baron et al. (2010) R. Baron, P. Boucaud, J. Carbonell, A. Deuzeman, V. Drach, et al., JHEP 1006, 111 (2010), eprint 1004.5284.
- Boucaud et al. (2008) P. Boucaud et al. (ETM), Comput. Phys. Commun. 179, 695 (2008), eprint 0803.0224.
- Bali et al. (2013) G. Bali, P. Bruns, S. Collins, M. Deka, B. Glasle, et al., Nucl.Phys. B866, 1 (2013), eprint 1206.7034.
- Alexandrou et al. (2013b) C. Alexandrou, M. Constantinou, V. Drach, K. Jansen, C. Kallidonis, et al. (2013b), eprint 1312.2874.
- Lüscher and Wolff (1990) M. Lüscher and U. Wolff, Nucl.Phys. B339, 222 (1990).