Extraction of multiple exponential signals from lattice correlation functions

We present a fast and simple algorithm that allows the extraction of multiple exponential signals from the temporal dependence of correlation functions evaluated on the lattice including the statistical fluctuations of each signal and treating properly backward signals. The basic steps of the method are the inversion of appropriate matrices and the determination of the roots of an appropriate polynomial, constructed using discretized derivatives of the correlation function. The method is tested strictly using fake data generated assuming a fixed number of exponential signals included in the correlation function with a controlled numerical precision and within given statistical fluctuations. All the exponential signals together with their statistical uncertainties are determined exactly by the algorithm. The only limiting factor is the numerical rounding off. In the case of correlation functions evaluated by large-scale QCD simulations on the lattice various sources of noise, other than the numerical rounding, can affect the correlation function and they represent the crucial factor limiting the number of exponential signals, related to the hadronic spectral decomposition of the correlation function, that can be properly extracted. The algorithm can be applied to a large variety of correlation functions typically encountered in QCD or QCD+QED simulations on the lattice, including the case of exponential signals corresponding to poles with arbitrary multiplicity and/or the case of oscillating signals. The method is able to to detect the specific structure of the multiple exponential signals without any a priori assumption and it determines accurately the ground-state signal without the need that the lattice temporal extension is large enough to allow the ground-state signal to be isolated.


S. Simula
Istituto Nazionale di Fisica Nucleare, Sezione di Roma Tre, Via della Vasca Navale 84, I-00146 Rome, Italy We present a fast and simple algorithm that allows the extraction of multiple exponential signals from the temporal dependence of correlation functions evaluated on the lattice including the statistical fluctuations of each signal and treating properly backward signals. The method starts from well-known features of the solution of ordinary (linear) differential equations (ODEs) and extracts multiple exponential signals from a generic correlation function simply by inverting appropriate matrices and by finding the roots of an appropriate polynomial, constructed using discretized derivatives of the correlation function. The method is tested strictly using fake data generated assuming a fixed number of exponential signals included in the correlation function with a controlled numerical precision and within given statistical fluctuations. All the exponential signals together with their statistical uncertainties are determined exactly by the ODE algorithm. The only limiting factor is the numerical rounding off. We show that, even when the total number of exponential signals contained in the correlation function is not known, the ODE method guarantees a quite good convergence toward accurate results for both masses and amplitudes, including their statistical fluctuations, at least for a significant subset of the exponential signals present in the correlation function. In the case of correlation functions evaluated by large-scale QCD simulations on the lattice various sources of noise, other than the numerical rounding, can affect the correlation function and they represent the crucial factor limiting the number of exponential signals, related to the hadronic spectral decomposition of the correlation function, that can be properly extracted. The ODE algorithm can be applied to a large variety of correlation functions typically encountered in QCD or QCD+QED simulations on the lattice, including the case of exponential signals corresponding to poles with arbitrary multiplicity and/or the case of oscillating signals. Among the appealing features of the ODE algorithm we mention its ability to detect the specific structure of the multiple exponential signals without any a priori assumption and the possibility to determine accurately the ground-state signal without the need that the lattice temporal extension is large enough to allow the ground-state signal to be isolated. the above sites are referred to as the source and the sink. Generally speaking the correlation function is integrated over the spatial extension of the lattice to get its dependence on the time distance t between the source and the sink. Since lattice simulations are defined in the Euclidean space, the temporal dependence of a correlation function admits a spectral decomposition in terms of a sum of exponential signals of the form Ae −M t , where M is a hadron mass (or energy when the hadron is moving), corresponding to an eigenvalue of the QCD (or QCD+QED) Hamiltonian, and A is the related amplitude containing a hadronic matrix element, which can be of interest.
In this work we present a fast and simple algorithm that allows the extraction of multiple exponential signals from the temporal dependence of correlation functions evaluated on the lattice including the determination of the statistical fluctuations of each signal. The method starts from well-known features of the solution of ordinary (linear) differential equations (ODEs) and extracts multiple exponential signals from a generic correlation function simply by inverting appropriate matrices and by finding the roots of an appropriate polynomial, constructed using discretized derivatives of the correlation function.
The idea of using properties of ODEs for extracting multiple exponential signals from the temporal dependence of lattice correlation functions is not at all a new one. Up to our knowledge variants of the method originally developed by Gaspard Riche de Prony long time ago [2], have been used to determine the masses of excited states from lattice correlators by finding the roots of an appropriate polynomial (see Refs. [3][4][5] and more recently Ref. [6]). Our method is characterized by the use of discretized derivatives of the correlation function.
The paper is organized as follows. In Section II the basic ingredients of the ODE algorithm are presented, namely the mass and amplitude matrices. The first step is the inversion of the mass matrix, which allows to construct the appropriate polynomial whose roots provide the masses of the exponential signals (i.e. the coefficients in the exponent). Subsequently, the amplitude matrix is constructed and inverted to determine the amplitudes of the exponential signals.
It will be shown that the ODE algorithm can be applied to a large variety of correlation functions typically encountered in large-scale QCD (or QCD+QED) simulations on the lattice, including the case of exponential signals corresponding to poles with arbitrary multiplicity and/or the case of oscillating signals. The ODE algorithm is able to detect the specific structure of the exponential signals (single/multiple poles or oscillating signals) without any a priori knowledge. This feature is a remarkable one, since it allows to specify the fitting ansatz appropriate for the lattice correlator.
The case of correlators with definite time parity is discussed in Section II A, while those corresponding to oscillating signals, poles with arbitrary multiplicity and multiple correlators are addressed in Sections II B, II C and II D, respectively. In Section II E we briefly discuss the issue of removing unwanted exponential signals from a correlator. Thanks to the ODE algorithm it is possible to filter out such signals without the need of determining their amplitudes, as described in details in the Appendix A.
In Section III the main features of the mass and amplitude matrices are discussed. Such matrices may be close to singularity and, therefore, the use of an appropriate numerical precision for matrix inversion is crucial for obtaining reliable results. The closeness to singularity can be described by means of the matrix condition number. On one hand side high values of the condition number represent a positive feature, that allows the ODE algorithm to be sensitive to the statistical fluctuations of the exponential signals. On the other hand side they may be a limiting factor related to the presence of noise in the correlation function, which can be produced by the numerical rounding and/or by other sources.
In Section III A the ODE algorithm is tested strictly using fake data for the correlation function, which are generated assuming a fixed number of exponential signals included in the correlation function and a controlled numerical precision. Each exponential signal is allowed to fluctuate within given uncertainties. We point out that, even if fake correlators represent an ideal situation different from the case of lattice correlators, it is nevertheless mandatory to check that the ODE algorithm (as well as any other algorithm developed for fitting the temporal dependence of lattice correlators) is able to provide exact results in a controlled situation. This is indeed our case and we show that the ODE algorithm is able to extract exactly all the exponential signals used as input together with their statistical uncertainties. An important, general feature of the ODE algorithm is that the ground-state signal can be extracted with accuracy even if the lattice temporal extension is not large enough to allow the ground-state signal to be isolated. This is a very useful property, which in particular can take care properly of the contamination of the excited states in the lattice correlators used for the determination of several hadronic quantities (like, e.g., the form factors).
In Section III B the case of correlators containing more exponential signals than those searched for is discussed. It will be shown how the ODE algorithm is still able to detect properly several exponential signals.
In Section III C the opposite situation, in which more signals than those included in the fake data are searched for, is investigated. It is found that the noise produced by the numerical rounding generates extra signals, which are not (or only weakly) suppressed exponentially in the time distance.
In Section IV the ODE algorithm is applied to the case of real data, namely correlation functions evaluated by means of large-scale QCD simulations on the lattice. Various sources of noise, other than the numerical rounding, can now affect the correlator, like e.g. the residues coming from gauge variant terms. It will be shown that the noise becomes the crucial factor limiting the number of exponentials, related to the eigenstates of the QCD Hamiltonian, that can be extracted from the correlation function.
In Section V the combination of the ODE algorithm with other techniques suitable for fitting lattice correlators is briefly discussed. An interesting possibility is represented by the combination of the ODE algorithm with a subsequent nonlinear least-squares minimizer, where masses and amplitudes are used as free parameters to be varied (without any prior) starting from the values obtained by the ODE method. In the case of the lattice correlators analyzed in Section IV the ODE solution is nicely confirmed, within the uncertainties, by the subsequent χ 2 -minimization procedure.
Finally, Section VI contains our conclusions.

II. THE ODE METHOD
Let's start by considering a correlator C(t) composed by N (+) exponential signals in the forward time direction and N (−) exponentials in the backward one: where T is the temporal extension of the lattice. In Eq.
(1) the masses M The correlator C(t) is supposed to be known at discretized values of the time distance t, namely t ≡ na with n = 1, ... N T , where a is the lattice spacing and N T ≡ T /a is the number of lattice points in the temporal direction. In lattice QCD (or QCD+QED) simulations a correlator of the form (1) may correspond, e.g., to the case of a nucleonic correlator, where the backward signals correspond to negative parity partners of the nucleon and its excitations.
Let's now consider the discretized (symmetric) time derivative 1 The generalization to the case of complex amplitudes is straightforward by adopting a strategy similar to the one described in Section II D.
where z m ≡ −sinh(a M m ) .
By repeating the application of the differential operation (5) one gets We also assume that the correlator C Note that, because of the presence of the factor (z m ) k in Eq. (7), at a fixed value of n the impact of the signals with higher masses increase as the order k of the derivative C (k) n increases.
The central step in our procedure consists in introducing N +1 real coefficients x k (k = 0, 1, ... N ) and considering the quantity where the polynomial P N (z) of degree N is given by The above polynomial has in general N roots depending on the coefficients x k . If the latter ones (which, we stress, are independent of n) are chosen so that the polynomial P N (z) has its roots at z = z m given by Eq. (6), then the condition holds for any value of n. Note that the roots z m of the polynomial P N (z) depend only on the masses M m and are independent of the amplitudes A m . Moreover, from Eq. (6) it follows that the roots z m are real numbers, positive for backward signals and negative for forward ones.
Equation (10) is a typical ordinary (linear) differential equation (ODE). Usually the coefficients x k are given and, therefore, the solution of Eq. (10) corresponds to Eq. (2) with the masses M m given by the roots (6) of the polynomial P N (z) and with the amplitudes A m depending on a suitable number of initial conditions.
Here we are interested in the inverse problem: starting from the known values of the correlator n and its derivatives we want to determine the coefficients x k of the polynomial (9) having its roots at z = z m . The procedure, which hereafter will be referred to as the ODE algorithm, is as follows.
Without any loss of generality we can put x N = 1 so that and Eq. (10) can be rewritten as The problem is to solve Eq. (12) for the N unknowns x k (k = 0, 1, ..., N − 1).
Our aim is to extract the multiple exponential signals in the correlator (2) using as input the into account. We get the following system of inhomogeneous linear equations where the N × N mass matrix M is given by and the vector V with dimension N by where σ (0) n is the uncertainty of the correlator (2), or where σ (k ) n is the uncertainty of the derivative (7). A more sophisticated choice is where D nn is the inverse of the covariance matrix of the correlator (2). We have checked that the performance of the ODE algorithm is not changed by the three choices (16-18) (see later Section III). In what follows we will use the definition (16). We point out that any autocorrelation between different values of n is taken into account by the use of the jackknife (bootstrap) procedure.
Thus, we rewrite Eqs. (14) and (15) as and Equations (19) and (20) are evaluated starting from the correlator (2) and its derivatives (7) for each jackknife or bootstrap event. Note that using the definition (16) for R (k ) n the system of linear equations (13) corresponds to minimize the variable χ 2 M defined as i.e. to the constraints ∂χ 2 M /∂x k = 0 for k = 0, 1, ...(N − 1) with x N = 1. For a non-singular matrix M the coefficients x k can be determined by inverting the matrix M : Once the coefficients x k are known, the roots of the polynomial P N (z) can be calculated, and thus the masses of the forward and backward exponential signals in lattice units, aM The last step is the determination of the amplitudes A m . To this end we introduce a χ 2 -variable defined as where σ (k) n is the statistical error of the derivative C (k) n . Then, we impose the minimization condition ∂χ 2 /∂ A m = 0, which leads to the following linear system of equations where The solution of the linear equation (25) is given by which allows to extract the forward and backward amplitudes, A Later in Sections III and IV we will check explicitly that the changes in the calculated amplitudes due to different choices of the χ 2 -variable (23) are totally negligible.
The key feature of the ODE method is the inversion of the mass and amplitude matrices, given respectively by Eqs. (19) and (26). This important issue will be discussed in Section III. In the following subsections we want to illustrate how the ODE method can be applied to specific forms of the correlation functions typically encountered in QCD (or QCD+QED) simulations on the lattice.
Let's consider the case in which the backward signals in the correlator (1) have the same masses and amplitudes (in absolute value) of those appearing in the forward signals, namely where (−) p i = ±1 is the parity with respect to the substitution t → (T − t). In what follows we will refer to (−) p i as the t-parity of the i-th exponential signal.
We can easily construct two combinations which have a definite t-parity: Thus, without any loss of generality we can consider a correlator C(t) of the following form where now (−) p can be equal to either −1 or +1 for all signals and N stands for the number of independent exponential signals, including now both the forward and the backward parts.
According to the results of the previous Section, the construction of the matrices (19) and (26) would imply working with matrices of dimension 2N × 2N . Instead, in the case of correlators with definite t-parity we can work with matrices having dimension N × N . To do that we modify the definition of the derivatives C (k) n . We now consider only even discretized derivatives and introduce the second derivative of the correlator (31) as where By repeating the application of the differential operation (32) one gets 2 2 We remind that the values of the correlator C signals. We write again explicitly the definitions of the mass and amplitudes matrices (and related vectors), since in the next Sections they will be repeatedly mentioned.
The constraint (10) is now replaced by where the coefficients x k for k = 0, 1, ..., (N − 1) can be determined by solving the linear system of with the N × N mass matrix given by and the N -dimensional vector V defined as By means of the coefficients x k we can construct the polynomial which has its roots at z = z i given by Eq. (33). Note that the roots z i : i) are positive for real values of M i , ii) do not depend on the amplitudes A i , and iii) are independent of the specific t-parity of the correlator (i.e., on the value of (−) p ). Note also that, thanks to the use of even where with σ (2k) n being the statistical error of the (even) derivative (34) and Note that the above functions depend on the specific t-parity of the correlator (31).
Before closing this subsection, we remind that at large time distances the signal of the lightest mass (e.g. M 1 ) is expected to dominate. Correspondingly the ratio C n may exhibit a plateau related to the value of aM 1 , namely We can therefore define an effective mass M

B. Oscillating signals
Specific formulations of the QCD action on the lattice may lead to quite specific features of the correlation functions constructed using quark and gluon propagators. One of such cases is represented by the staggered formulation, in which the correlation functions may have oscillating terms related to the presence of opposite (spatial) parity partners.
In order to simplify the notations we limit ourselves to the case of a correlator with positive t-parity, namely where now N (±) stands for normal/oscillating exponential signals.
The possible presence of oscillating signals do not represent a problem for the ODE algorithm.
We can therefore apply our ODE algorithm. The only difference is that now the roots z m , while remaining real, can be either positive or negative. The positive ones correspond to normal signals, while the negative roots to oscillating signals (more precisely z m ≤ −4 because Before closing this subsection, we point out that the possibility to detect the presence of oscillating signals occurs when the analyzed correlator has a definite time parity, i.e. when the relation between masses and roots is given by Eq. (33), which originates from the use of even derivatives only.
On the contrary, when both odd and even derivatives are used, the relation between masses and roots is given by Eq. (6). Since sinh(aM i ± iπ) = −sinh(aM i ) = sinh(−aM i ), a forward(backward) oscillating signal cannot be distinguished from a backward(forward) non-oscillating signal, basing only on the analyses of the mass matrix (14) constructed using both odd and even derivatives.

C. Poles with arbitrary multiplicity
Till now we have considered multiple exponential signals corresponding to single poles of the In this subsection we address the case of poles characterized by an arbitrary multiplicity µ i , which correspond to exponential signals multiplied by a polynomial of degree (µ i − 1) in the time distance t. Thus, the correlator (31) is replaced by where N ≡ N i=1 µ i represents the total number of exponential terms [i.e., those in the square brackets in the r.h.s. of Eq. (52)]. In modern lattice QCD+QED simulations the above situation may occur, e.g., when isospin breaking effects, due to the quark electric charges and to the mass difference δm between u and d quarks, are taken into account at leading order in the electromagnetic coupling α em and in δm (see, e.g., Ref. [7]). In this case correlation functions contain double poles (i.e. µ i = 2).
Within the ODE algorithm the procedure is as follows. Let's start from the correlator where the amplitudes B iµ ≡ B iµ a µ are given in lattice units and The sequence of even derivatives can be constructed for k = 1, ..., N . One gets x k z k and z i is given by Eq. (33). The constraint which means that the root z i of the polynomial P N (z) has multiplicity µ i , namely (with x N = 1)  In this subsection we address briefly the case in which different correlators sharing the same masses M i are available. In lattice QCD (or QCD+QED) simulations such a situation may occur when different interpolating fields, like e.g. local (L) and smeared (S) fields, are adopted in the source and in the sink.
For sake of simplicity let's consider the simple case of four correlators C f f (t) with f (f ) = L, S.
We look for multiple exponential signals of the form We assume also that for each jackknife (or bootstrap) event one has C LS (t) = C SL (t) (if not, one can average over the two types of correlators or use the one with the smallest statistical fluctuations). Thus in what follows we limit ourselves to consider three independent correlators: C LL (t), C LS (t) and C SS (t) sharing the same masses M i and, respectively, with real amplitudes where the coefficients x k should be independent of the specific choice of the correlator C f f (t).
In order to guarantee such an independence and to make a simultaneous use of the information contained in all the correlators we construct the variable χ 2 M given by and we impose the global constraints ∂ χ 2 M /∂x k = 0 for k = 0, 1, ...(N − 1) with x N = 1. It follows that the coefficients x k should satisfy the linear system of equations where consequently of the common masses M i . Finally, the amplitudes A L i and A S i can be easily obtained after applying the procedure illustrated in Section II A to the individual correlators C f f (t).

E. Filtering
Lattice correlation functions are defined in a (discretized) Euclidean space after performing the Wick rotation from the physical Minkowsky space. After the rotation, however, correlation functions may contain exponential signals with masses lighter than the one relevant for the physical process under investigation. A well-known example is the K → ππ decay, where according to Ref. [8] the energy non-conserving matrix elements with the state of two pions at rest are involved.
In these cases, since at finite lattice volume the eigenvalues of the QCD Hamiltonian are discretized, a possible procedure is to try to subtract or, in other words, to filter out the unwanted signals.
The filtering of unwanted signals can be carried out after the determination of both masses and amplitudes for all the exponential signals present in the correlator. In Appendix A we describe a simple procedure, based on the ODE algorithm, that allows to filter out exponential signals from a given correlator without the need of determining the amplitudes of the unwanted signals.

III. INVERSION OF THE MASS AND AMPLITUDE MATRICES
The numerical inversion of the mass and amplitude matrices can be carried out using standard methods, like the lower-upper decomposition (or factorization) method [9]. In what follows we will refer to the mass and amplitude matrices as defined in Section II A.
If the interval of summation over n in Eq. (37) is restricted to a single value of n, the matrix M has a vanishing determinant and therefore it cannot be inverted. Generally speaking, the sum over various values of n protects against the singularity of the matrix M . However, the mass matrix may still be close to singularity. As we shall see in this Section and in Section IV, the quasi-singularity represents on one hand side a positive feature, that allows the ODE algorithm to be sensitive to the fluctuations of the multiple exponential signals, and on the other hand side a limiting factor related to the presence of noise in the correlation function.
In the field of numerical analysis a condition number κ(M ) can be associated to the system of linear equations (36) and provides a bound on the relative error of the solution (35) with respect to the relative error of the vector (38). The condition number κ(M ) is a property of the matrix M and it is defined as [10] where the symbol ||M || stands for the norm of the matrix M . The latter can be defined in several ways and hereafter we adopt the Frobenius definition [9] ||M || ≡ If the matrix M is singular, then κ(M ) → ∞. It is easy to show that where δx (δV ) is the error on x (V ) and the vector norm is consistently defined as ||x|| =  [11], which allows to change the precision level during run time. We have used at least 32 digits of precision (quadruple precision) reaching 96 digits in the most severe cases.
Also the amplitude matrix (41) needs to be inverted. In this case the condition number κ(A) turns out to be much smaller and it can be handled properly by adopting 32 digits of precision.
For finding the roots of the polynomial P N (z) we make use of the powerful, open-source software package MPSolve [12].

A. Fake data
We now test the ODE algorithm by generating fake data for the correlator C (0) n to be used as benchmarks. We stress that, even if fake correlators represent an ideal situation different from the case of lattice correlators, it is nevertheless mandatory to check that our algorithm (as well as any other algorithm developed for fitting the temporal dependence of lattice correlators) is able to provide exact results in a controlled situation.
Let us start by considering a correlator with positive t-parity of the form on a lattice with temporal extension N T = T /a = 96 and spacing a (the specific value of a is not required since we work in lattice units). The fake data are generated by allowing the amplitudes A i and the masses aM i to fluctuate with uncertainties δA i and δ(aM i ), respectively, adopting (uncorrelated) Gaussian distributions to produce a total of 40 jackknives.
The number of exponential signals in Eq. (68) is taken to be equal to N = 12. The values chosen for the amplitudes and the masses are collected in Table I   In the case at hand, for the relative deviations we get We have also considered the cases in which the uncertainties δ(aM i ) and δA i are either increased up to 10% or decreased down to 0.1% of the corresponding masses aM i and amplitudes A i . In the former case the maximum relative deviations are ∆ max ∼ 2 · 10 −13 and δ∆ max ∼ 3 · 10 −11 , while in the latter case we get ∆ max ∼ 3 · 10 −13 and δ∆ max ∼ 2 · 10 −9 . can be detected, since they lead to large variations of the function f (x).
As more exponential signals with masses above the heaviest one in Table I are added in the fake correlator (68), both the condition numbers κ(M ) and κ(A) as well as the maximum relative deviations ∆ max and δ∆ max increase quickly as shown in Table II.
12 ∼ 2 · 10 −13 ∼ 2 · 10 −10 ≈ 10 31 ∼ 3 · 10 5 14 ∼ 2 · 10 −11 ∼ 4 · 10 −8 ≈ 10 38 ∼ 5 · 10 6 16 ∼ 2 · 10 −4 ∼ 6 · 10 −2 ≈ 10 49 ∼ 1 · 10 9  The quick rise of ∆ max and δ∆ max with N is related to the impact of the numerical rounding of the fake correlator, while it does not depend on the numerical precision of the ODE method, which can be always kept at the desired level. Indeed, the numerical accuracy of the fake correlator C (0) n is external to the ODE method and it governs the maximum number of exponential signals that can be determined precisely (i.e. within given values of the maximum relative deviations ∆ max and δ∆ max ). In Table III we have collected the values of N max found indicatively for different levels of the numerical precision of the correlator C (0) n while keeping ∆ max < 5 · 10 −9 and δ∆ max < 5 · 10 −6 .
We now want to discuss a particular set of masses and amplitudes for the correlator (68), which is relevant for simulating a typical situation occurring when hadronic form factors are extracted from appropriate lattice correlators. The usual procedure is to form a suitable ratio of correlation functions (typically, the ratio of 3-point and 2-point correlation functions), which at large time distances exhibits a plateau. From the latter the hadronic form factor is obtained. We stress that the above procedure requires that the temporal separation between the source and the sink should be large enough to allow the ground-state signal to be isolated. The contamination of excited states is usually investigated by varying the separation, which however may be computationally quite expensive.
In order to simulate the above situation we put the lightest mass in the fake correlator (68) equal to zero, obtaining in this way a correlator that becomes constant at large time distances. In  Fig. 2, indicates clearly that the dominance of the ground-state signal, having a vanishing mass, is not yet reached. This is not at all a problem for the ODE algorithm, which is able to determine accurately all the masses and amplitudes of Table IV. We point out that the above results illustrate an important, general feature of the ODE algorithm: it allows to extract accurately the ground-state signal without the need that the lattice temporal extension (or the temporal separation between the source and the sink for 3-point correlators) is large enough to allow the ground-state signal to be isolated. In other words, the ODE method is able to remove properly the contamination of excited states also at relatively small values of the lattice temporal extension (or without changing the temporal separation between the source and the sink for 3-point correlators).
Let us now consider the case of a correlator containing poles with arbitrary multiplicity, namely where µ i is the multiplicity of the i-th exponential with i = 1, 2, ...N and N i=1 µ i = N . The values chosen for the masses aM i , the amplitudes A iµ and the multiplicities µ i are collected in Table V together with the corresponding uncertainties δ(aM i ) and δA iµ . There are four single poles, two double poles and a quadruple pole for a total of N = 12 exponential signals. The relative uncertainties are taken to be equal to 1% in the case of the masses and 5% for the amplitudes.
The quadruple precision (32 digits) is used for evaluating the fake correlator (71) and the octuple precision (64 digits) for the internal ODE calculations. We now apply the ODE algorithm using all the derivatives C An interesting case is represented by the following fake correlator composed by six double poles with masses and amplitudes given in Table VI Table V. imaginary part equal to −π). The numerical precision for generating the fake correlator (72) is again 32 digits (quadruple precision). The condition numbers of the mass and amplitude matrices are ∼ 5·10 30 and ∼ 3·10 10 , respectively.
The accuracy of the ODE results for all the masses, amplitudes and their uncertainties is better than ∼ 1 ppb.
Results with the same quality can be obtained in the case of multiple correlators (see Section II D), which we do not report here for sake of brevity. We just mention that a well-established procedure to deal with multiple correlators is represented by the method based on the Generalized Eigenvalue Problem (GEVP) [13]. With this method it is possible to extract masses and amplitudes of both ground and excited states. In particular, with four correlators, generated using two different interpolating fields, the ground and the first excited states can be determined. As the number of excited states increases, the number of interpolating fields and correspondingly the number of correlators should be increased. This is at variance with the ODE method, which is able to detect properly many exponential signals independently on the number of correlators used.

B. ODE analyses with N ODE < N
In this Section we address the case in which the total number of derivatives C (2k) n included in the construction of the mass and amplitude matrices is less than the total number of exponential signals included in the fake data for the correlator C with masses and amplitudes given in Table I for N = 12.
We apply the ODE algorithm using C with k, k = 0, 1, ..., (N ODE − 1) and the vector (38) is now given by The solution of the ODE conditions provides the coefficients x k , which are used to construct the polynomial of degree N ODE where Generally speaking we expect that the ODE masses aM ODE j represent an approximation of the lighter N ODE masses aM j included in the correlator C This is indeed the phenomenology we observe. In Tables VII and VIII we have collected the results obtained by applying the ODE algorithm assuming N ODE = 8 in the ranges n = [5,92] (n 0 = 4) and n = [10, 87] (n 0 = 9), respectively. First of all it can be seen that in Table 5 Table I  in the fake data of the correlator C Though not all the ODE masses and amplitudes reproduce the input values, the fake data of the correlator C (0) n are reasonably reproduced by its ODE representation We find that the absolute value of the residue [C The number of exponential signals properly reproduced by the ODE algorithm increases as N ODE increases. As already pointed out in Section II, this is due to the fact that in any fixed range of values of n the derivative C (2k) n is more sensitive to the signals with higher masses as the order k increases (cf., e.g., the factor (z i ) k in Eq. (34)). This feature is illustrated in Fig. 6,  Table I), versus the value of n 0 + 1, which defines the range of the analysis n = [n 0 + 1, N T − n 0 ].
The results for the various states are slightly shifted horizontally for better readability. multiplicity µ j greater than 1. One has with s j n given in Eq. (82) and (−) p being the t-parity of the correlator.
Finally, in the ODE procedure we introduce two tolerance factors δ R and δ I , which help in finding multiple poles. They are defined as follows. When the real parts of the masses of two signals are close enough, namely the two signals are replaced by a signal having as real part [µ j aRe(M ODE j )+µ j aRe(M ODE j )]/(µ j + µ j ) and multiplicity µ j + µ j . When the imaginary part of the mass of a signal is small enough, it is put to zero. By using the above tolerance factors the relative error in the ODE representation (84)) is of the order O(δ 2 ).
The results presented in this Section illustrate that, even when the total number of exponential signals contained in the fake data is not known, the ODE method guarantees a quite good convergence toward accurate results for both masses and amplitudes, including their statistical fluctuations, at least for a significant subset of the exponential signals present in the fake correlator. We stress that the specific structure of the ODE representation [see Eq. (83) or Eq. (84)] does not require any a priori assumption, but it is properly detected by the ODE method. According to the numerical precision of the correlator the ODE algorithm detects almost exactly all the masses and the amplitudes (together with their statistical fluctuations) used as input (see Table III).
In this Section we describe an interesting feature of the ODE algorithm when N ODE > N , namely when the number of roots of the polynomial P N ODE (z) (see Eq. (76)) is larger than the number N of signals present in the correlator C n .
Let us consider again the fake data for the correlator (68) corresponding only to the first eight masses and amplitudes given in Table I   trying to take into account the small residual terms. Note that the extra masses of Table IX exhibit an approximate pattern, whose interpretation requires however a separate study. The noisy roots can be simply discarded while keeping the remaining eight roots. The ODE amplitudes A ODE j corresponding to the non-noisy signals can be determined from the same amplitude matrix calculated when N ODE = N = 8.

IV. ANALYSIS OF LATTICE CORRELATORS
In this Section we apply the ODE algorithm to selected correlation functions evaluated by means of large-scale QCD simulations on the lattice. We will employ in particular correlators evaluated using gauge ensembles produced by the European (now Extended) Twisted Mass Collaboration (ETMC). The numerical precision of the available correlators is optimistically the double precision (16 digits). From Table III we expect that the ODE algorithm should detect accurately no more than 6 exponential signals. However, in the case of lattice correlators the noise is not only due to the numerical rounding, but e.g. also to the residues coming from gauge variant terms. Therefore, we expect that only 3-4 exponential signals can be identified properly by the ODE algorithm.
Let's start by considering the 2-point correlator C P S (t) defined as where P 5 (x) = q 2 (x)γ 5 q 1 (x) is a local interpolating field that creates in x a pseudoscalar (PS) mesons made of two valence quarksq 1 and q 2 with masses m 1 and m 2 . At large time distances one has so that the meson mass M P S and the matrix element Z P S = P S|q 2 γ 5 q 1 |0 can be extracted from the exponential fit given in the r.h.s. of Eq. (89).
The evaluation of the correlator (88) involves the determination of the so-called all-to-all quark propagator. For the latter the statistical accuracy is significantly improved by using the so-called "one-end" stochastic method [14], which includes spatial stochastic sources at a single time slice chosen randomly.
Among the gauge ensembles generated by ETMC with N f = 2 + 1 + 1 dynamical quarks (see In what follows we refer to the non-noisy states as the physical states. where N phys indicates the number of physical states. The difference C P S (t) − C phys (t) is therefore related to the noisy states only and it is shown in the right panel of Fig. 7 for N ODE = 6, 8, 10.  Table X. A good convergence for the ground state and progressively for the excited states is observed as N ODE increases. For N ODE ≥ 6 a nice agreement for both masses and amplitudes occurs within the uncertainties. The latter ones appear to be slightly larger as N ODE increases above N ODE = 6 and this is probably due to the occurrence of more extra roots in the mass matrix.
Note the strong stability of the mass and amplitude for the ground-state signal for N ODE ≥ 6.
The approach commonly used to extract the ground-state signal is to look at the time dependence of the effective mass in order to identify the range of values of n, where the ground-state dominates. In Fig. 8 the results corresponding to two definitions of the effective mass are shown.
One definition is given by Eq. (45) and the other one is the usual logarithmic slope   We now move to the heavier sector of charmonium. We select the gauge ensemble D20.48 [15] corresponding to a lattice volume V × T = 48 3 × 96a 4 and to a lattice spacing a 0.062 fm, and choose the valence quark masses equal to m 1 = m 2 1.18 GeV (in the MS(2 GeV) scheme), i.e. close to the physical charm quark mass [15]. The number of independent gauge configurations employed is 100. We consider both the PS correlator (88) and the vector (V) one given by where V i (x) = q 2 (x)γ i q 1 (x) is a local interpolating field that creates a vector meson in x. Four stochastic sources per gauge configuration are employed leading to a relative statistical error which does not exceed ∼ 0.4% and ∼ 0.7% for the PS and V correlators, respectively.
We apply the ODE algorithm with N ODE ≥ 4 finding always four physical signals and N ODE −4 noisy states (all with multiplicity equal to 1). The results for the physical signals obtained with N ODE = 4, 6 and 8 are collected in Tables XI and XII for the PS and V correlators, respectively.
The convergence of both masses and amplitudes for all the four physical states is quite good.
The time dependence of the effective masses (45) and (91) is shown in Fig. 9 Note that the ODE method does not correspond to the minimization of a unique χ 2 -variable.
Indeed, as pointed out in Section II, the inversion of the mass matrix is equivalent to minimize the variable χ 2 M defined by Eq. (21), while, once the masses are given, the determination of the amplitudes corresponds to the minimization of the χ 2 -variable given by Eqs. (23-24) (or just its first term with k = 0). Furthermore, as shown in the previous Section, the ODE algorithm is sensitive to the presence of noisy states only through the mass matrix and independently of the size of the corresponding amplitudes. Instead, a χ 2 -minimization procedure is expected to be sensitive also to the amplitudes of the noise. In this way the combined ODE plus χ 2 -minimization procedure can produce either a non-trivial check of the ODE solution or possibly a refinement.
We have explicitly used the above combination in the case of the lattice correlators analyzed in Section IV. We have found that the ODE solution is nicely confirmed, within the uncertainties, by the subsequent χ 2 -minimization procedure. An interesting feature is that no priors on the masses and amplitudes are required in the χ 2 minimization in order to obtain stable results. In the constrained curve fitting method of Ref. [17] priors are instead introduced just for stabilizing the results of the fitting procedure.
The combination of the ODE method with nonlinear least-squares minimizers requires, however, dedicated numerical investigations, which are outside the scope of the present paper.
We close this Section by observing that the number of hadronic states in a lattice QCD correlator (i.e. the number of exponential signals) is not finite and, generally speaking, it increases as the time distance decreases. In this respect it is worth to mention the representation of the vector-vector current correlator, relevant for the determination of the hadronic contribution to the anomalous magnetic moment of the muon, developed using quark-hadron duality at short time distances in Ref. [16]. The dual contribution represents an effective way to perform a resummation of an infinite number of highly excited hadronic states at short time distances. Thus, the ODE algorithm can be combined with quark-hadron duality: it can be applied not to the full correlator, but to its difference with the dual contribution. Such an issue will be investigated in a separate work.

VI. CONCLUSIONS
We have presented a fast and simple algorithm that allows the extraction of multiple exponential signals from the temporal dependence of correlation functions evaluated on the lattice including the statistical fluctuations of each signal and treating properly backward signals.
The method starts from well-known features of the solution of ordinary (linear) differential equations and extracts multiple exponential signals from a generic correlation function simply by inverting appropriate mass and amplitude matrices and by finding the roots of an appropriate polynomial. The method is based on the use of discretized derivatives of the correlation function.
An important feature of the ODE method is the level of singularity of the mass matrix, described by its condition number. On one hand side this represents a positive feature, that allows the ODE algorithm to be sensitive to the fluctuations of the exponential signals, and on the other hand side it is a limiting factor related to the presence of noise in the correlation function. Huge values of the condition number require an appropriate treatment of the numerical precision, which has been obtained adopting the multiple precision software from Ref. [11]. Moreover, the root finding has been carried out accurately using the open-source software package MPSolve [12].
We have tested extensively the ODE method using fake data, generated assuming a fixed number of exponential signals included in the correlator with a controlled numerical precision and within given statistical fluctuations. All the exponential signals with their statistical uncertainties are determined exactly by the ODE algorithm, when the total number of exponential signals is known.
The only limiting factor is the numerical rounding off. We have shown that, even when the total number of exponential signals contained in the correlator is not known, the ODE method guarantees a quite good convergence toward accurate results for both masses and amplitudes, including their statistical fluctuations, at least for a significant subset of the exponential signals present in the correlator.
Then, few cases of correlation functions evaluated by means of large-scale QCD simulations on the lattice have been addressed explicitly. In the case of lattice correlators, several sources of noise, other than the numerical rounding, can affect the correlator. As shown in Section IV, the noise represents the crucial factor limiting the number of physical exponential signals, related to the hadronic spectral decomposition of the correlation function, that can be determined.
We have illustrated that the ODE algorithm can be applied to a large variety of correlation functions typically encountered in QCD or QCD+QED simulations on the lattice, including the case of exponential signals corresponding to poles with arbitrary multiplicity and/or the case of oscillating signals. Two important features of the ODE algorithm are: i) its ability to detect the proper structure of the multiple exponential signals without any a priori assumption, and ii) the extraction of the ground-state signal with accuracy without the need that the lattice temporal extension is large enough to allow the ground-state signal to be isolated. This is a very useful property, which in particular can take care properly of the contamination of the excited states in the lattice correlators used for the determination of hadronic quantities, like e.g. the form factors.
A further application of the ODE algorithm is represented by its combination with a subsequent nonlinear least-squares minimizer, where masses and amplitudes are used as free parameters (without any prior) to be varied starting from the values obtained by the ODE method.
A careful study of the origin and the structure of the noisy states, which are systematically In this Appendix we make use of the ODE algorithm to address the issue of subtracting exponential signals from a given correlator without the need of determining their amplitudes.
Without loss of generality let's consider the case of a correlator with a given t-parity (−) p , composed by N exponential signals each with multiplicity equal to 1, namely A i e −aM i n + (−) p e −aM i (N T −n) .
The first step of the ODE algorithm is to calculate iteratively the derivatives with We now apply the ODE algorithm to the modified correlator C The generalization to the case of roots with arbitrary multiplicities is straightforward.
We point out that the above filtering procedure can be very easily adapted when the masses of the unwanted signals are known a priori or from other analyses not based on the ODE algorithm. In such cases the coefficients y k of the polynomial Q M (z) can be calculated directly without applying the ODE algorithm to the original correlator (A1). Then, the filtering of the unwanted states can proceed by constructing the modified correlator (A9) and by applying to it the ODE algorithm.