Extraction of multiple exponential signals from lattice correlation functions

TL;DR

This paper introduces a fast algorithm for extracting multiple exponential signals from lattice correlation functions, effectively handling statistical fluctuations and backward signals, with applications in QCD simulations.

Contribution

The paper presents a novel, efficient method for extracting multiple exponential signals from lattice correlation functions without prior assumptions, improving analysis of spectral data.

Findings

Successfully tested with fake data, accurately extracting signals and uncertainties.

Limited primarily by numerical rounding and noise in real lattice data.

Applicable to various correlation functions in QCD and QCD+QED simulations.

Abstract

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.

Tables

Table 1. Table 1: Values of the masses a ​ M i 𝑎 subscript 𝑀 𝑖 aM_{i} , amplitudes A i subscript 𝐴 𝑖 A_{i} and their uncertainties δ ​ ( a ​ M i ) 𝛿 𝑎 subscript 𝑀 𝑖 \delta(aM_{i}) and δ ​ A i 𝛿 subscript 𝐴 𝑖 \delta A_{i} adopted for generating the fake data for the correlator ( 68 ) assuming N = 12 𝑁 12 N=12 . The temporal extension is N T = T / a = 96 subscript 𝑁 𝑇 𝑇 𝑎 96 N_{T}=T/a=96 .

$i$	$a{M}_i$	$A_i$	$i$	$a{M}_i$	$A_i$
$1$	$0.05\pm 0.0005$	$0.30\pm 0.0030$	$7$	$0.95\pm 0.0095$	$0.30\pm 0.0030$
$2$	$0.01\pm 0.0010$	$0.70\pm 0.0070$	$8$	$1.10\pm 0.0110$	$-0.45\pm 0.0045$
$3$	$0.25\pm 0.0025$	$-0.10\pm 0.0010$	$9$	$1.30\pm 0.0130$	$0.50\pm 0.0050$
$4$	$0.40\pm 0.0040$	$0.80\pm 0.0080$	$10$	$1.55\pm 0.0155$	$0.30\pm 0.0030$
$5$	$0.65\pm 0.0065$	$-0.30\pm 0.0030$	$11$	$1.80\pm 0.0180$	$0.20\pm 0.0020$
$6$	$0.75\pm 0.0075$	$0.25\pm 0.0025$	$12$	$2.10\pm 0.0210$	$-0.55\pm 0.0055$

Table 2. Table 2: Values of the the maximum relative deviations Δ m ​ a ​ x subscript Δ 𝑚 𝑎 𝑥 \Delta_{max} and δ ​ Δ m ​ a ​ x 𝛿 subscript Δ 𝑚 𝑎 𝑥 \delta\Delta_{max} versus the number N 𝑁 N of exponential signals included in the fake correlator C n ( 0 ) superscript subscript 𝐶 𝑛 0 C_{n}^{(0)} . The numerical precision of C n ( 0 ) superscript subscript 𝐶 𝑛 0 C_{n}^{(0)} corresponds to 32 32 32 digits (quadruple precision), while the internal numerical precision adopted in the ODE algorithm corresponds to 64 64 64 digits (octuple precision). The temporal extension of the lattice is N T = T / a = 96 subscript 𝑁 𝑇 𝑇 𝑎 96 N_{T}=T/a=96 and the uncertainties δ ​ ( a ​ M i ) 𝛿 𝑎 subscript 𝑀 𝑖 \delta(aM_{i}) and δ ​ A i 𝛿 subscript 𝐴 𝑖 \delta A_{i} are equal to 1 % percent 1 1\% of the corresponding masses and amplitudes a ​ M i 𝑎 subscript 𝑀 𝑖 aM_{i} and A i subscript 𝐴 𝑖 A_{i} .

$N$	${\Delta }_{max}$	$\delta {\Delta }_{max}$	$\kappa \left(\overline{M}\right)$	$\kappa \left(\overline{A}\right)$
$12$	$\sim 2\cdot {10}^{-13}$	$\sim 2\cdot {10}^{-10}$	$\approx {10}^{31}$	$\sim 3\cdot {10}^5$
$14$	$\sim 2\cdot {10}^{-11}$	$\sim 4\cdot {10}^{-8}$	$\approx {10}^{38}$	$\sim 5\cdot {10}^6$
$16$	$\sim 2\cdot {10}^{-4}$	$\sim 6\cdot {10}^{-2}$	$\approx {10}^{49}$	$\sim 1\cdot {10}^9$

Table 3. Table 3: Values of the the maximum number N m ​ a ​ x subscript 𝑁 𝑚 𝑎 𝑥 N_{max} of exponential signals that can be determined by the ODE algorithm keeping Δ m ​ a ​ x < 5 ⋅ 10 − 9 subscript Δ 𝑚 𝑎 𝑥 ⋅ 5 superscript 10 9 \Delta_{max}<5\cdot 10^{-9} and δ ​ Δ m ​ a ​ x < 5 ⋅ 10 − 6 𝛿 subscript Δ 𝑚 𝑎 𝑥 ⋅ 5 superscript 10 6 \delta\Delta_{max}<5\cdot 10^{-6} for different levels of the numerical precision of the correlator C n ( 0 ) superscript subscript 𝐶 𝑛 0 C_{n}^{(0)} . The third, fourth and fifth columns contain the condition numbers κ ​ ( M ¯ ) 𝜅 ¯ 𝑀 \kappa(\overline{M}) , κ ​ ( A ¯ ) 𝜅 ¯ 𝐴 \kappa(\overline{A}) and the internal numerical precision adopted in the ODE algorithm, respectively. The temporal extension of the lattice is N T = T / a = 96 subscript 𝑁 𝑇 𝑇 𝑎 96 N_{T}=T/a=96 and the uncertainties δ ​ ( a ​ M i ) 𝛿 𝑎 subscript 𝑀 𝑖 \delta(aM_{i}) and δ ​ A i 𝛿 subscript 𝐴 𝑖 \delta A_{i} are equal to 1 % percent 1 1\% of the corresponding masses and amplitudes a ​ M i 𝑎 subscript 𝑀 𝑖 aM_{i} and A i subscript 𝐴 𝑖 A_{i} .

correlator precision	$N_{max}$	$\kappa \left(\overline{M}\right)$	$\kappa \left(\overline{A}\right)$	internal precision
(digits)				(digits)
$16$	$6$	$\approx {10}^{14}$	$\sim 3\cdot {10}^2$	$32$
$32$	$14$	$\approx {10}^{38}$	$\sim 5\cdot {10}^6$	$64$
$64$	$20$	$\approx {10}^{74}$	$\sim 7\cdot {10}^{14}$	$96$

Table 4. Table 4: Values of the masses a ​ M i 𝑎 subscript 𝑀 𝑖 aM_{i} , amplitudes A i subscript 𝐴 𝑖 A_{i} and their uncertainties δ ​ ( a ​ M i ) 𝛿 𝑎 subscript 𝑀 𝑖 \delta(aM_{i}) and δ ​ A i 𝛿 subscript 𝐴 𝑖 \delta A_{i} adopted for generating the fake data for the correlator ( 68 ) assuming N = 4 𝑁 4 N=4 . The temporal extension is N T = T / a = 48 subscript 𝑁 𝑇 𝑇 𝑎 48 N_{T}=T/a=48 .

$i$	$a{M}_i$	$A_i$
$1$	$0.000\pm 0.000$	$0.500\pm 0.025$
$2$	$0.100\pm 0.001$	$0.800\pm 0.040$
$3$	$0.200\pm 0.002$	$1.000\pm 0.050$
$4$	$0.500\pm 0.005$	$1.500\pm 0.075$

Table 5. Table 5: Values of the masses a ​ M i 𝑎 subscript 𝑀 𝑖 aM_{i} , amplitudes A i ​ μ subscript 𝐴 𝑖 𝜇 A_{i\mu} , multiplicities μ i subscript 𝜇 𝑖 \mu_{i} and the uncertainties δ ​ ( a ​ M i ) 𝛿 𝑎 subscript 𝑀 𝑖 \delta(aM_{i}) and δ ​ A i ​ μ 𝛿 subscript 𝐴 𝑖 𝜇 \delta A_{i\mu} adopted for generating the fake data for the correlator ( 71 ). The temporal extension is N T = T / a = 96 subscript 𝑁 𝑇 𝑇 𝑎 96 N_{T}=T/a=96 .

$i$	$a{M}_i$	${\mu }_i$	$A_{i\mu }$	$i$	$a{M}_i$	${\mu }_i$	$A_{i\mu }$
$1$	$\mathrm{\hspace{0.17em}\hspace{0.17em}0.15}\pm 0.0015$	$2$	$\mathrm{\hspace{0.17em}\hspace{0.17em}0.050}\pm 0.0025$	$5$	$1.20\pm 0.012$	$4$	$\mathrm{\hspace{0.17em}\hspace{0.17em}0.50}\pm 0.025$
			$\mathrm{\hspace{0.17em}\hspace{0.17em}0.002}\pm 0.0010$				$\mathrm{\hspace{0.17em}\hspace{0.17em}0.60}\pm 0.030$
$2$	$0.30\pm 0.003$	$1$	$-0.10\pm 0.005$				$-0.20\pm 0.010$
$3$	$0.60\pm 0.006$	$2$	$-0.20\pm 0.010$				$\mathrm{\hspace{0.17em}\hspace{0.17em}1.00}\pm 0.050$
			$\mathrm{\hspace{0.17em}\hspace{0.17em}0.500}\pm 0.025$	$6$	$1.50\pm 0.015$	$1$	$\mathrm{\hspace{0.17em}\hspace{0.17em}1.50}\pm 0.075$
$4$	$0.90\pm 0.009$	$1$	$\mathrm{\hspace{0.17em}\hspace{0.17em}1.000}\pm 0.050$	$7$	$2.00\pm 0.020$	$1$	$\mathrm{\hspace{0.17em}\hspace{0.17em}2.00}\pm 0.100$

Table 6. Table 6: Values of the masses a ​ M ~ m 𝑎 subscript ~ 𝑀 𝑚 a\widetilde{M}_{m} and the amplitudes A ~ m subscript ~ 𝐴 𝑚 \widetilde{A}_{m} and B ~ m subscript ~ 𝐵 𝑚 \widetilde{B}_{m} ( m = 1 , … , 6 𝑚 1 … 6 m=1,...,6 ) together with their uncertainties adopted for generating the fake data for the correlator ( 72 ). The temporal extension is N T = T / a = 96 subscript 𝑁 𝑇 𝑇 𝑎 96 N_{T}=T/a=96 .

$m$	$a{\tilde{M}}_m$	${\tilde{A}}_m$	${\tilde{B}}_m$
$1$	$0.05\pm 0.0020$	$1.00\pm 0.0020$	$0.50\pm 0.0040$
$2$	$0.20\pm 0.0030-i\pi$	$0.70\pm 0.0030$	$-0.20\pm 0.0060$
$3$	$0.55\pm 0.0020$	$0.80\pm 0.0030$	$-0.30\pm 0.0060$
$4$	$0.80\pm 0.0040-i\pi$	$0.25\pm 0.0080$	$0.30\pm 0.0060$
$5$	$1.20\pm 0.0050$	$-0.45\pm 0.0100$	$0.50\pm 0.0120$
$6$	$1.60\pm 0.0020-i\pi$	$0.30\pm 0.0140$	$0.30\pm 0.0140$

Table 7. Table 7: Values of the masses a ​ M j O ​ D ​ E 𝑎 superscript subscript 𝑀 𝑗 𝑂 𝐷 𝐸 aM_{j}^{ODE} , amplitudes A j O ​ D ​ E superscript subscript 𝐴 𝑗 𝑂 𝐷 𝐸 A_{j}^{ODE} and their uncertainties obtained by applying the ODE algorithm with N O ​ D ​ E = 8 subscript 𝑁 𝑂 𝐷 𝐸 8 N_{ODE}=8 to the correlator ( 68 ) containing N = 12 𝑁 12 N=12 exponential signals with masses and amplitudes given in Table  1 . The range of the ODE analysis is n = [ n 0 + 1 , N T − n 0 ] = [ 5 ,  92 ] 𝑛 subscript 𝑛 0 1 subscript 𝑁 𝑇 subscript 𝑛 0 5 92 n=[n_{0}+1,\,N_{T}-n_{0}]=[5,\,92] . All the roots of the polynomial ( 76 ) are found to have multiplicity equal to μ i = 1 subscript 𝜇 𝑖 1 \mu_{i}=1

$j$	$a{M}_j^{ODE}$	$A_j^{ODE}$
$1$	$0.04995\pm 0.00050$	$0.29883\pm 0.00307$
$2$	$0.09986\pm 0.00100$	$0.70001\pm 0.00697$
$3$	$0.25781\pm 0.00323$	$-0.11321\pm 0.00421$
$4$	$0.39552\pm 0.00418$	$0.79056\pm 0.00911$
$5$	$0.72065\pm 0.00522+i\left(0.110\pm 0.025\right)$	$0.05472\pm 0.01392$
$6$	$0.72065\pm 0.00522-i\left(0.110\pm 0.025\right)$	$-0.02925\pm 0.00472$
$7$	$1.5015\pm 0.01899$	$0.67175\pm 0.03405$
$8$	$2.1090\pm 0.05427$	$-0.43229\pm 0.01820$

Table 8. Table 8: The same as in Table  7 , but for a range of the ODE analysis equal to n = [ n 0 + 1 , N T − n 0 ] = [ 10 ,  87 ] 𝑛 subscript 𝑛 0 1 subscript 𝑁 𝑇 subscript 𝑛 0 10 87 n=[n_{0}+1,\,N_{T}-n_{0}]=[10,\,87] . In the fourth and fifth columns the values of the masses and amplitudes for the first eight signals of Table  1 are shown for ease of comparison.

$j$	$a{M}_j^{ODE}$	$A_j^{ODE}$	$a{M}_j$	$A_j$
$1$	$0.05000\pm 0.00050$	$0.30000\pm 0.00300$	$0.05\pm 0.0005$	$0.30\pm 0.0030$
$2$	$0.10000\pm 0.00100$	$0.70000\pm 0.00700$	$0.01\pm 0.0010$	$0.70\pm 0.0070$
$3$	$0.24993\pm 0.00250$	$-0.09987\pm 0.00100$	$0.25\pm 0.0025$	$-0.10\pm 0.0010$
$4$	$0.40009\pm 0.00401$	$0.80080\pm 0.00800$	$0.40\pm 0.0040$	$0.80\pm 0.0080$
$5$	$0.63611\pm 0.00575$	$-0.22323\pm 0.01413$	$0.65\pm 0.0065$	$-0.30\pm 0.0030$
$6$	$0.82329\pm 0.00451$	$0.30990\pm 0.02311$	$0.75\pm 0.0075$	$0.25\pm 0.0025$
$7$	$1.17580\pm 0.05828$	$-0.21524\pm 0.08739$	$0.95\pm 0.0095$	$0.30\pm 0.0030$
$8$	$1.40830\pm 0.01651$	$0.72549\pm 0.08468$	$1.10\pm 0.0110$	$-0.45\pm 0.0045$

Table 9. Table 9: Masses a ​ M j O ​ D ​ E 𝑎 superscript subscript 𝑀 𝑗 𝑂 𝐷 𝐸 aM_{j}^{ODE} of the extra ( N O ​ D ​ E − 8 subscript 𝑁 𝑂 𝐷 𝐸 8 N_{ODE}-8 ) signals found by the ODE algorithm when N O ​ D ​ E = 10 , 12 , 14 , 16 subscript 𝑁 𝑂 𝐷 𝐸 10 12 14 16 N_{ODE}=10,12,14,16 is used for analyzing the correlator ( 68 ), in which eight exponential signals, corresponding to the first eight masses and amplitudes given in Table  1 , have been included.

$j$	$N_{ODE}=10$	$N_{ODE}=12$	$N_{ODE}=14$	$N_{ODE}=16$
$9$	$i\left(2.873\pm 0.495\right)$	$i\left(2.960\pm 0.477\right)$	$i\left(2.996\pm 0.446\right)$	$i\left(3.023\pm 0.479\right)$
$10$	$i\left(2.312\pm 0.561\right)$	$i\left(2.576\pm 0.428\right)$	$i\left(2.722\pm 0.472\right)$	$i\left(2.804\pm 0.433\right)$
$11$		$i\left(2.211\pm 0.431\right)$	$i\left(2.430\pm 0.472\right)$	$i\left(2.535\pm 0.439\right)$
$12$		$i\left(1.799\pm 0.486\right)$	$i\left(2.138\pm 0.432\right)$	$i\left(2.317\pm 0.487\right)$
$13$			$i\left(1.822\pm 0.375\right)$	$i\left(2.089\pm 0.437\right)$
$14$			$i\left(1.482\pm 0.409\right)$	$i\left(1.804\pm 0.413\right)$
$15$				$i\left(1.554\pm 0.430\right)$
$16$				$i\left(1.273\pm 0.460\right)$

Table 10. Table 10: Masses a ​ M j O ​ D ​ E 𝑎 superscript subscript 𝑀 𝑗 𝑂 𝐷 𝐸 aM_{j}^{ODE} (upper panel) and amplitudes A j O ​ D ​ E superscript subscript 𝐴 𝑗 𝑂 𝐷 𝐸 A_{j}^{ODE} (lower panel) of the four physical signals found by the ODE algorithm when N O ​ D ​ E = 4 subscript 𝑁 𝑂 𝐷 𝐸 4 N_{ODE}=4 , 5 5 5 , 6 6 6 , 8 8 8 and 10 10 10 is used for analyzing the PS correlator ( 88 ) calculated in the case of the ETMC gauge ensemble B25.32  Carrasco:2014cwa for valence quark masses equal to m 1 = m 2 = m ℓ ≃ 12 subscript 𝑚 1 subscript 𝑚 2 subscript 𝑚 ℓ similar-to-or-equals 12 m_{1}=m_{2}=m_{\ell}\simeq 12 MeV in the MS ¯ ​ ( 2 ​ GeV ) ¯ MS 2 GeV \overline{\rm MS}(2~{}\mbox{GeV}) scheme.

$j$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$
	$\left(N_{ODE}=4\right)$	$\left(N_{ODE}=5\right)$	$\left(N_{ODE}=6\right)$	$\left(N_{ODE}=8\right)$	$\left(N_{ODE}=10\right)$
$1$	$0.10935\pm 0.00088$	$0.10644\pm 0.00075$	$0.10693\pm 0.00027$	$0.10694\pm 0.00034$	$0.10697\pm 0.00034$
$2$	$\mathrm{\hspace{0.17em}0.9625}\pm 0.0452$	$\mathrm{\hspace{0.17em}0.5126}\pm 0.0903$	$\mathrm{\hspace{0.17em}0.5839}\pm 0.0307$	$\mathrm{\hspace{0.17em}0.5600}\pm 0.0463$	$\mathrm{\hspace{0.17em}0.5872}\pm 0.0432$
$3$	$\mathrm{\hspace{0.17em}2.0201}\pm 0.0311$	$\mathrm{\hspace{0.17em}1.1453}\pm 0.0931$	$\mathrm{\hspace{0.17em}1.2147}\pm 0.0586$	$\mathrm{\hspace{0.17em}1.1570}\pm 0.0920$	$\mathrm{\hspace{0.17em}1.2168}\pm 0.1004$
$4$		$\mathrm{\hspace{0.17em}2.0910}\pm 0.0568$	$\mathrm{\hspace{0.17em}2.1321}\pm 0.0427$	$\mathrm{\hspace{0.17em}2.0888}\pm 0.0651$	$\mathrm{\hspace{0.17em}2.1311}\pm 0.0819$

Table 11. Table 11: Masses a ​ M j O ​ D ​ E 𝑎 superscript subscript 𝑀 𝑗 𝑂 𝐷 𝐸 aM_{j}^{ODE} (upper panel) and amplitudes A j O ​ D ​ E superscript subscript 𝐴 𝑗 𝑂 𝐷 𝐸 A_{j}^{ODE} (lower panel) of the four physical signals found by the ODE algorithm when N O ​ D ​ E = 4 subscript 𝑁 𝑂 𝐷 𝐸 4 N_{ODE}=4 , 6 6 6 and 8 8 8 is used for analyzing the PS correlator ( 88 ) calculated in the case of the ETMC gauge ensemble D20.48  Carrasco:2014cwa for valence quark masses in the charm region, namely m 1 = m 2 ≃ 1.18 subscript 𝑚 1 subscript 𝑚 2 similar-to-or-equals 1.18 m_{1}=m_{2}\simeq 1.18 GeV in the MS ¯ ​ ( 2 ​ GeV ) ¯ MS 2 GeV \overline{\rm MS}(2~{}\mbox{GeV}) scheme.

$j$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$
	$\left(N_{ODE}=4\right)$	$\left(N_{ODE}=6\right)$	$\left(N_{ODE}=8\right)$
$1$	$0.95201\pm 0.00015$	$0.95192\pm 0.00015$	$0.95189\pm 0.00019$
$2$	$\mathrm{\hspace{0.17em}1.2466}\pm 0.0054$	$\mathrm{\hspace{0.17em}1.2348}\pm 0.0035$	$\mathrm{\hspace{0.17em}1.2304}\pm 0.0051$
$3$	$\mathrm{\hspace{0.17em}1.7938}\pm 0.0121$	$\mathrm{\hspace{0.17em}1.7589}\pm 0.0093$	$\mathrm{\hspace{0.17em}1.7480}\pm 0.0130$
$4$	$\mathrm{\hspace{0.17em}2.6731}\pm 0.0284$	$\mathrm{\hspace{0.17em}2.5969}\pm 0.0213$	$\mathrm{\hspace{0.17em}2.5784}\pm 0.0247$

Table 12. Table 12: The same as in Table  11 but in the case of the vector correlator ( 92 ).

$j$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$	$a{M}_j^{ODE}$
	$\left(N_{ODE}=4\right)$	$\left(N_{ODE}=6\right)$	$\left(N_{ODE}=8\right)$
$1$	$0.99431\pm 0.00104$	$0.99556\pm 0.00027$	$0.99550\pm 0.00022$
$2$	$\mathrm{\hspace{0.17em}1.2458}\pm 0.0243$	$\mathrm{\hspace{0.17em}1.2706}\pm 0.0070$	$\mathrm{\hspace{0.17em}1.2689}\pm 0.0058$
$3$	$\mathrm{\hspace{0.17em}1.8276}\pm 0.0330$	$\mathrm{\hspace{0.17em}1.8472}\pm 0.0153$	$\mathrm{\hspace{0.17em}1.8431}\pm 0.0134$
$4$	$\mathrm{\hspace{0.17em}2.7898}\pm 0.0354$	$\mathrm{\hspace{0.17em}2.7943}\pm 0.0217$	$\mathrm{\hspace{0.17em}2.7892}\pm 0.0197$

Table 13. Table 13: Ground-state mass, a ​ M 1 𝑎 subscript 𝑀 1 aM_{1} , and amplitude, A 1 subscript 𝐴 1 A_{1} , extracted from a single exponential fit in a chosen range n = [ n m ​ i ​ n , n m ​ a ​ x ] 𝑛 subscript 𝑛 𝑚 𝑖 𝑛 subscript 𝑛 𝑚 𝑎 𝑥 n=[n_{min},\,n_{max}] in the case of the V correlator ( 92 ) corresponding to two valence quarks with masses m 1 = m 2 ≃ 1.18 subscript 𝑚 1 subscript 𝑚 2 similar-to-or-equals 1.18 m_{1}=m_{2}\simeq 1.18 GeV in the case of the ETMC gauge ensemble D20.48  Carrasco:2014cwa .

$[n_{min},n_{max}]$	$a{M}_1$	$A_1$
$[25,\mathrm{\hspace{0.17em}45}]$	$0.99492\pm 0.00018$	$0.01642\pm 0.00013$
$[30,\mathrm{\hspace{0.17em}45}]$	$0.99463\pm 0.00024$	$0.01625\pm 0.00018$
$[25,\mathrm{\hspace{0.17em}40}]$	$0.99494\pm 0.00020$	$0.01643\pm 0.00013$
$[30,\mathrm{\hspace{0.17em}40}]$	$0.99460\pm 0.00022$	$0.01623\pm 0.00018$

Equations

$$C(t)=\sum _{i=1}^{N^{(+)}}{A}_i^{(+)}{e}^{-M_i^{(+)}t}+\sum _{j=1}^{N^{(-)}}{A}_j^{(-)}{e}^{-M_j^{(-)}(T-t)},$$

$$C(t=na)\equiv C_n^{(0)}=\sum _{m=1}^N{\tilde{A}}_m{e}^{-a{\tilde{M}}_{m}n}$$

$${\tilde{M}}_m=M_i^{(+)},{\tilde{A}}_m=A_i^{(+)}$$

$${\tilde{M}}_m=-M_j^{(-)},{\tilde{A}}_m=A_j^{(-)}{e}^{-M_j^{(-)}T}$$

$$C_n^{(1)}\equiv \frac{1}{2}\left[C_{n+1}^{(0)}-C_{n-1}^{(0)}\right]=\sum _{m=1}^N{\tilde{A}}_m{z}_m{e}^{-a{\tilde{M}}_{m}n},$$

$$z_m\equiv -\text{mbox}sinh(a{\tilde{M}}_m).$$

$$C_n^{(k)}=\frac{1}{2}\left[C_{n+1}^{(k-1)}-C_{n-1}^{(k-1)}\right]=\sum _{m=1}^N{\tilde{A}}_m(z_m{)}^k{e}^{-a{\tilde{M}}_{m}n}.$$

$$\sum _{k=0}^N{x}_k{C}_n^{(k)}=\sum _{m=1}^N{\tilde{A}}_m\left[\sum _{k=0}^N{x}_k{z}_m^k\right]e^{-a{\tilde{M}}_{m}n}=\sum _{m=1}^N{\tilde{A}}_m{P}_N(z_m)e^{-a{\tilde{M}}_{m}n},$$

$$P_N(z)\equiv \sum _{k=0}^N{x}_k{z}^k.$$

$$\sum _{k=0}^N{x}_k{C}_n^{(k)}=0$$

$$P_N(z)\equiv \sum _{k=0}^N{x}_k{z}^k=\sum _{k=0}^{N-1}{x}_k{z}^k+z^N=\prod _{m=1}^N(z-z_m)$$

$$\sum _{k=0}^{N-1}{x}_k{C}_n^{(k)}=-C_n^{(N)}.$$

$$\sum _{k=0}^{N-1}{M}_{k^{\prime }k}{x}_k=V_{k^{\prime }},$$

$$M_{k^{\prime }k}\equiv \sum _{n=N+1}^{N_T-N}{R}_n^{(k^{\prime })}{C}_n^{(k)},$$

$$V_{k^{\prime }}\equiv -\sum _{n=N+1}^{N_T-N}{R}_n^{(k^{\prime })}{C}_n^{(N)}.$$

$$R_n^{(k^{\prime })}=\frac{C_n^{(k^{\prime })}}{[{\sigma }_n^{(0)}{]}^2},$$

$$R_n^{(k^{\prime })}=\frac{C_n^{(k^{\prime })}}{[{\sigma }_n^{(k^{\prime })}{]}^2},$$

$$R_n^{(k^{\prime })}=D_{n{n}^{\prime }}^{(0)}{C}_{n^{\prime }}^{(k^{\prime })},$$

$$M_{k^{\prime }k}\equiv \sum _{n=N+1}^{N_T-N}\frac{C_n^{(k^{\prime })}{C}_n^{(k)}}{[{\sigma }_n^{(0)}{]}^2},$$

$$V_{k^{\prime }}\equiv -\sum _{n=N+1}^{N_T-N}\frac{C_n^{(k^{\prime })}{C}_n^{(N)}}{[{\sigma }_n^{(0)}{]}^2}.$$

$${\chi }_M^2\equiv \sum _{n=N+1}^{N_T-N}\frac{1}{[{\sigma }_n^{(0)}{]}^2}{\left[\sum _{k^{\prime }=0}^N{x}_{k\prime }{C}_n^{(k^{\prime })}\right]}^2,$$

$$x_k=\sum _{k^{\prime }=0}^{N-1}{M}_{k{k}^{\prime }}^{-1}{V}_{k^{\prime }}.$$

$${\chi }^2$$

$${\chi }_k^2$$

$$\sum _{m^{\prime }=1}^N{A}_{m{m}^{\prime }}{\tilde{A}}_{m^{\prime }}=W_m,$$

$$A_{m{m}^{\prime }}$$

$$W_m$$

$${\tilde{A}}_m=\sum _{m^{\prime }=1}^N{A}_{m{m}^{\prime }}^{-1}{W}_{m^{\prime }},$$

$$C(t)=\sum _{i=1}^N{A}_i\left[e^{-M_{i}t}+(-{)}^{p_i}{e}^{-M_i(T-t)}\right],$$

$$C^{(\pm )}(t)=\frac{1}{2}\left[C(t)\pm C(T-t)\right]=\sum _{i=1}^N{A}_i\left[e^{-M_{i}t}\pm e^{-M_i(T-t)}\right].$$

$$C(t)=\sum _{i=1}^N{A}_i\left[e^{-M_{i}t}+(-{)}^p{e}^{-M_i(T-t)}\right],$$

Full text

I Introduction

The investigation of the limits of the Standard Model of particle physics and the search for possible signatures of New Physics represent nowadays the main task for both experimental and theoretical physics. As present (and planned) experimental facilities explore new energy frontiers and improve the precision of the measurements, the importance of flavor physics is continuously growing. It is crucial for such studies to quantify the non-perturbative effects due to the strong interaction among hadrons in physical processes. Large-scale numerical simulations of Quantum Chromodynamics (QCD) performed on the lattice allow for an ab initio computation of hadronic quantities based on first principles only. Since several years various quantities relevant for flavour physics phenomenology have been determined by lattice QCD simulations reaching the impressive level of precision of ${\cal{O}}(1\%)$ or even better (see the recent FLAG-4 review Aoki:2019cca ). Moreover, in recent years lattice computations include also isospin-breaking effects due to electromagnetism and to the mass difference between up and down quarks, which are of the order of ${\cal{O}}(\alpha_{em})$ and ${\cal{O}}[(m_{d}-m_{u})/\Lambda_{\mathrm{QCD}}]$, i.e. both of the order of ${\cal{O}}(1\%)$. Thus, QCD+QED simulations on the lattice are now crucial for making further progresses in flavor physics phenomenology.

The main quantities evaluated through QCD (or QCD+QED) simulations on the lattice are correlation functions defined in terms of the expectation value on the vacuum of the time-ordered product of operators appropriate for the investigation of the physical process of interest. The operators are located at different sites on the lattice. In the case of 2-point correlation functions 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^{-Mt}$, 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 data is not at all a new one and it traces back to the method originally developed by Gaspard Riche de Prony long time ago Prony . Since then, there have been many developments and applications aimed at fitting data with linear combinations of real (or complex) exponentials within many disciplines in Science and Engineering (see, e.g., Ref. PS and references therein). In the last decades Prony-type methods PS ; OS have been used to determine the masses of excited states from lattice correlators (see Refs. Fleming:2004hs ; Fleming:2009wb ; Beane:2009kya and more recently Ref. Cushman:2019tcv ). These methods are based on the construction of a difference equation which involves data at equally spaced values in the time direction and the coefficients of the exponents are determined by finding the roots of a characteristic polynomial. Our method, which we have developed having in mind applications to the analysis of lattice correlation functions, is a Prony-type method characterized by the direct 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.1, while those corresponding to oscillating signals, poles with arbitrary multiplicity and multiple correlators are addressed in Sections II.2, II.3 and II.4, respectively. In Section II.5 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.1 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.2 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.3 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 $\chi^{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:

[TABLE]

where $T$ is the temporal extension of the lattice. In Eq. (1) the masses $M_{i}^{(+)}$ and $M_{j}^{(-)}$ are nonnegative real numbers (as in the case of hadronic masses) and the amplitudes $A_{i}^{(+)}$ and $A_{j}^{(-)}$ are real numbers111The generalization to the case of complex amplitudes is straightforward by adopting a strategy similar to the one described in Section II.4..

The correlator $C(t)$ is supposed to be known at discretized values of the time distance $t$, namely $t\equiv na$ with $n=1,...\,N_{T}$, where $a$ is the lattice spacing and $N_{T}\equiv 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.

For sake of simplicity we will refer to the quantities $M_{i}^{(+)}$ and $M_{j}^{(-)}$ as masses. It is however clear that Eq. (1) may correspond also to the case of correlation functions for moving hadrons by simply replacing hadron masses with energies.

The correlator (1) can be rewritten as

[TABLE]

with $N\equiv N^{(+)}+N^{(-)}$ and

[TABLE]

in the case of forward signals ($m=i=1,...\,N^{(+)}$) and

[TABLE]

for backward signals ($m=N^{(+)}+j=N^{(+)}+1,...\,N^{(+)}+N^{(-)}$).

Let’s now consider the discretized (symmetric) time derivative

[TABLE]

where

[TABLE]

By repeating the application of the differential operation (5) one gets

[TABLE]

We also assume that the correlator $C_{n}^{(0)}$ is known, for each time distance, in terms of a number of jackknife or bootstrap events. Its statistical error $\sigma_{n}^{(0)}$ can be correspondingly evaluated. Starting from the correlator $C_{n}^{(0)}$, the sequence of the correlators $C_{n}^{(k)}$ with $k=1,...\,N$ can be evaluated for each jackknife or bootstrap event together with their statistical errors $\sigma_{n}^{(k)}$. It is understood that what follows applies for each single jackknife or bootstrap event.

The values of the correlator $C_{n}^{(0)}$ are provided in the range $n=[1,\,N_{T}]$, while the derivatives $C_{n}^{(k)}$ for $k=1,...\,N$ are evaluated only in the range $n=[k+1,\,N_{T}-k]$. Outside this range the values of the derivatives are not independent and, therefore, we put $C_{n}^{(k)}=0$.

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_{n}^{(k)}$ increases.

The central step in our procedure consists in introducing $N+1$ real coefficients $x_{k}$ ($k=0,1,...\,N$) and considering the quantity

[TABLE]

where the polynomial $P_{N}(z)$ of degree N is given by

[TABLE]

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

[TABLE]

holds for any value of $n$. Note that the roots $z_{m}$ of the polynomial $P_{N}(z)$ depend only on the masses $\widetilde{M}_{m}$ and are independent of the amplitudes $\widetilde{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 $\widetilde{M}_{m}$ given by the roots (6) of the polynomial $P_{N}(z)$ and with the amplitudes $\widetilde{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 $C_{n}^{(0)}$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

[TABLE]

and Eq. (10) can be rewritten as

[TABLE]

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 knowledge of the correlator in a given range of values of $n$, which eventually can span the full temporal extension $[1,\,N_{T}]$. Therefore, we multiply Eq. (12) by a set of $N$ functions $R_{n}^{(k^{\prime})}$ with $k^{\prime}=0,1,...\,(N-1)$ and sum over $n$ in a given range from $n_{min}$ to $n_{max}$. Since the largest range in which the derivatives $C_{n}^{(k)}$ can be calculated is $n=[k+1,\,N_{T}-k]$, we put directly $n_{min}=N+1$ and $n_{max}=N_{T}-N$, so that all the values of the correlator (2) in the full range $[1,\,N_{T}]$ are taken into account. We get the following system of inhomogeneous linear equations

[TABLE]

where the $N\times N$ mass matrix $M$ is given by

[TABLE]

and the vector $V$ with dimension $N$ by

[TABLE]

The choice of the functions $R_{n}^{(k^{\prime})}$ is in principle arbitrary, provided it leads to a non-singular mass matrix $M$. We have explored different choices for $R_{n}^{(k^{\prime})}$. Simple and natural choices are either

[TABLE]

where $\sigma_{n}^{(0)}$ is the uncertainty of the correlator (2), or

[TABLE]

where $\sigma_{n}^{(k^{\prime})}$ is the uncertainty of the derivative (7). A more sophisticated choice is

[TABLE]

where $D_{nn^{\prime}}^{(0)}$ 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

[TABLE]

and

[TABLE]

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_{n}^{(k^{\prime})}$ the system of linear equations (13) corresponds to minimize the variable $\chi_{M}^{2}$ defined as

[TABLE]

i.e. to the constraints $\partial\chi_{M}^{2}/\partial 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$:

[TABLE]

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_{i}^{(+)}$ and $aM_{j}^{(-)}$, can be determined from Eq. (6).

The last step is the determination of the amplitudes $\widetilde{A}_{m}$. To this end we introduce a $\chi^{2}$-variable defined as

[TABLE]

where $\sigma_{n}^{(k)}$ is the statistical error of the derivative $C_{n}^{(k)}$. Then, we impose the minimization condition $\partial\chi^{2}/\partial\widetilde{A}_{m}=0$, which leads to the following linear system of equations

[TABLE]

where

[TABLE]

The solution of the linear equation (25) is given by

[TABLE]

which allows to extract the forward and backward amplitudes, $A_{i}^{(+)}$ and $A_{j}^{(-)}$ using Eqs. (3-4).

We stress that the choice of the $\chi^{2}$-variable given by Eqs. (23-24) is not unique and, consequently, also the definitions of the amplitude matrix (26) and vector (27). For instance, one can limit the sum in Eq. (23) to the first term $k=0$ without involving the derivatives $C_{n}^{(k)}$ with $k>0$. Correspondingly also in Eqs. (26-27) the sum over $k$ should be limited to the first term $k=0$ only. Later in Sections III and IV we will check explicitly that the changes in the calculated amplitudes due to different choices of the $\chi^{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.

II.1 Correlators with definite time parity

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

[TABLE]

where $(-)^{p_{i}}=\pm 1$ is the parity with respect to the substitution $t\to(T-t)$. In what follows we will refer to $(-)^{p_{i}}$ as the $t$-parity of the $i$-th exponential signal. Note that, because of periodicity, each exponential signal in Eq. (29) receives contributions also from multiple wrappings around the lattice time extension. The latter ones correspond to a modification of the amplitudes $A_{i}$ by a multiplicative factor $(1-e^{-M_{i}T})^{-1}$, which we consider already absorbed into the amplitude $A_{i}$ appearing in Eq. (29). The multiplicative factor is numerically small for enough large values of $M_{i}T$, as it happens in many applications of lattice QCD simulations.

We can easily construct two combinations which have a definite $t$-parity:

[TABLE]

Thus, without any loss of generality we can consider a correlator $C(t)$ of the following form

[TABLE]

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\times 2N$. Instead, in the case of correlators with definite $t$-parity we can work with matrices having dimension $N\times N$. To do that we modify the definition of the derivatives $C_{n}^{(k)}$. We now consider only even discretized derivatives and introduce the second derivative of the correlator (31) as

[TABLE]

where

[TABLE]

By repeating the application of the differential operation (32) one gets222We remind that the values of the correlator $C_{n}^{(0)}$ are provided in the range $n=[1,\,N_{T}]$, while the derivatives $C_{n}^{(2k)}$ for $k=1,..\,N$ can be determined only in the range $n=[k+1,\,N_{T}-k]$. Outside this range we consider that $C_{n}^{(2k)}=0$.

[TABLE]

We now repeat the same steps followed previously in the case of different forward and backward 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

[TABLE]

where the coefficients $\overline{x}_{k}$ for $k=0,1,...,(N-1)$ can be determined by solving the linear system of equations (choosing $\overline{x}_{N}=1$)

[TABLE]

with the $N\times N$ mass matrix given by

[TABLE]

and the $N$-dimensional vector $\overline{V}$ defined as

[TABLE]

By means of the coefficients $\overline{x}_{k}$ we can construct the polynomial

[TABLE]

which has its roots at $z=\overline{z}_{i}$ given by Eq. (33). Note that the roots $\overline{z}_{i}$:   i) are positive for real values of $M_{i}$,   ii) do not depend on the amplitudes $A_{i}$ (i.e., the determination of the nonlinear unknowns $M_{i}$ is independent on the values of the linear ones $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 derivatives only [see the r.h.s. of Eq. (34)], the constraint (35) can be satisfied simultaneously by the forward and backward parts of each signal.

Following the same procedure described in the previous Section [see Eqs. (23-27)] the amplitudes $A_{i}$ can be obtained by solving the linear system of equations

[TABLE]

where

[TABLE]

with $\sigma_{n}^{(2k)}$ being the statistical error of the (even) derivative (34) and

[TABLE]

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}^{(2)}/C_{n}^{(0)}$ may exhibit a plateau related to the value of $aM_{1}$, namely

[TABLE]

We can therefore define an effective mass $M_{eff}^{(cosh)}$ as

[TABLE]

II.2 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

[TABLE]

where now $N^{(\pm)}$ stands for normal/oscillating exponential signals.

The possible presence of oscillating signals do not represent a problem for the ODE algorithm. Indeed, since $(-)^{n}=e^{\pm i\pi n}$ the correlator (46) can be written as the sum of $N=N^{(+)}+N^{(-)}$ exponential signals as

[TABLE]

where

[TABLE]

in the case of the normal signals ($m=i=1,...\,N^{(+)}$) and

[TABLE]

for oscillating signals ($m=N^{(+)}+j=N^{(+)}+1,...\,N^{(+)}+N^{(-)}$), for which the masses $a\widetilde{M}_{m}$ acquire an imaginary part equal to $-\pi$ (or equivalently $+\pi$).

We can therefore apply our ODE algorithm. The only difference is that now the roots $\overline{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 $\overline{z}_{m}\leq-4$ because $2\left[\mbox{cosh}(aM-i\pi)-1\right]=-2\left[\mbox{cosh}(aM)+1\right]$).

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 $\mbox{sinh}(aM_{i}\pm i\pi)=-\mbox{sinh}(aM_{i})=\mbox{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.

II.3 Poles with arbitrary multiplicity

Till now we have considered multiple exponential signals corresponding to single poles of the form $(p^{2}-M_{i}^{2})^{-1}$ in (Minkowskian) momentum space. In this subsection we address the case of poles characterized by an arbitrary multiplicity $\mu_{i}$, which correspond to exponential signals multiplied by a polynomial of degree ($\mu_{i}-1$) in the time distance $t$. Thus, the correlator (31) is replaced by

[TABLE]

where $N\equiv\sum_{i=1}^{\overline{N}}\mu_{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 $\delta m$ between $u$ and $d$ quarks, are taken into account at leading order in the electromagnetic coupling $\alpha_{em}$ and in $\delta m$ (see, e.g., Ref. RM123 ). In this case correlation functions contain double poles (i.e. $\mu_{i}=2$).

Within the ODE algorithm the procedure is as follows. Let’s start from the correlator

[TABLE]

where the amplitudes $\overline{B}_{i\mu}\equiv B_{i\mu}a^{\mu}$ are given in lattice units and

[TABLE]

The sequence of even derivatives

[TABLE]

can be constructed for $k=1,...,N$. One gets

[TABLE]

where $\overline{P}_{N}(z)\equiv\sum_{k=0}^{N}\overline{x}_{k}z^{k}$ and $\overline{z}_{i}$ is given by Eq. (33). The constraint $\sum_{k=0}^{N}\overline{x}_{k}C_{n}^{(2k)}=0$ is satisfied for each value of $n$ only if

[TABLE]

which means that the root $\overline{z}_{i}$ of the polynomial $\overline{P}_{N}(z)$ has multiplicity $\mu_{i}$, namely (with $\overline{x}_{N}=1$)

[TABLE]

Thus, the masses $M_{i}$ and their multiplicities $\mu_{i}$ ($i=1,2,...,M$) are determined in two steps: first by solving the linear system of equations corresponding to Eqs. (36-38) to obtain the coefficients $x_{k}$ and then by finding the roots of the polynomial $\overline{P}_{N}(z)\equiv\sum_{k=0}^{N-1}\overline{x}_{k}z^{k}+z^{N}$ with their multiplicities $\mu_{i}$. Finally, the amplitudes $\overline{B}_{i\mu}$ can be determined by solving the linear system of equations corresponding to Eqs. (40-42) with the functions $f_{n}^{(i)}$ replaced by $f_{n}^{(i,\mu-0)}$ given in Eq. (54).

II.4 Multiple correlators

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^{ff^{\prime}}(t)$ with $f(f^{\prime})=L,S$. We look for multiple exponential signals of the form

[TABLE]

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 $A_{i}^{L}A_{i}^{L}$, $A_{i}^{L}A_{i}^{S}$ and $A_{i}^{S}A_{i}^{S}$.

The ODE algorithm can be applied separately to each correlator $C^{ff^{\prime}}(t)$, namely Eq. (36) becomes

[TABLE]

where the coefficients $\overline{x}_{k}$ should be independent of the specific choice of the correlator $C^{ff^{\prime}}(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 $\widehat{\chi}_{M}^{2}$ given by

[TABLE]

and we impose the global constraints $\partial\widehat{\chi}_{M}^{2}/\partial\overline{x}_{k}=0$ for $k=0,1,...(N-1)$ with $\overline{x}_{N}=1$. It follows that the coefficients $\overline{x}_{k}$ should satisfy the linear system of equations

[TABLE]

where

[TABLE]

The solution of Eq. (62) allows the determination of the roots of the polynomial (39) and consequently of the common masses $M_{i}$ . Finally, the amplitudes $A_{i}^{L}$ and $A_{i}^{S}$ can be easily obtained after applying the procedure illustrated in Section II.1 to the individual correlators $C^{ff^{\prime}}(t)$.

II.5 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\to\pi\pi$ decay, where according to Ref. Maiani:1990ca 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 LU . In what follows we will refer to the mass and amplitude matrices as defined in Section II.1.

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 $\kappa(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 $\kappa(M)$ is a property of the matrix $M$ and it is defined as cond

[TABLE]

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 LU

[TABLE]

If the matrix $M$ is singular, then $\kappa(M)\to\infty$. It is easy to show that

[TABLE]

where $\delta x$ ($\delta V$) is the error on $x$ ($V$) and the vector norm is consistently defined as $||x||=\sqrt{\sum_{i=1}^{N}x_{i}^{2}}$. Thus, as a rule of thumb, when the condition number $\kappa(M)$ is equal to $\approx 10^{d}$, one may lose up to $d$ digits of accuracy in solving numerically Eq. (36).

The condition number $\kappa(M)$ depends strongly on the dimension $N$ of the matrix $M$. We have found that the appropriate handling of the numerical inversion of Eq. (36) requires the use of multiple precision software. In this work we have adopted the open-source software package MPFUN2015 mpfun , 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 $\kappa(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 MPSolve .

III.1 Fake data

We now test the ODE algorithm by generating fake data for the correlator $C_{n}^{(0)}$ 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

[TABLE]

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 $\delta A_{i}$ and $\delta(aM_{i})$, respectively, adopting (uncorrelated) Gaussian distributions to produce a total of 40 jackknives. We remind that the ODE procedure is always applied to each single jackknife. The corresponding results provide central value and errors for the masses and the amplitudes according to the jackknife procedure.

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 1 together with their uncertainties $\delta A_{i}$ and $\delta(aM_{i})$. The latter ones are taken to be equal to $1\%$ of the corresponding amplitudes $A_{i}$ and masses $aM_{i}$.

The time dependencies of the correlator $C_{n}^{(0)}$, of the derivatives $C_{n}^{(2k)}$ for $k=1,2,3$ and of the effective mass $aM_{eff}^{(cosh)}$ (see Eq. (45)) are shown in Fig. 1. The values of the lightest few masses have been chosen in such a way that at large time distances the ratio $C_{n}^{(2)}/C_{n}^{(0)}$ does not exhibit any plateaux within half of the temporal extension of the lattice. Thus, no direct information on the lightest mass can be extracted from the fake correlator using the standard effective mass methodology.

The numerical precision for generating the fake correlator (68) with $N=12$ is chosen to be 32 digits, i.e. quadruple precision (see later on Table 3), while all the internal ODE numerical calculations, namely the evaluation of the derivatives $C_{n}^{(2k)}$ for $k=1,...12$ as well as the numerical inversion of the mass and amplitude matrices of dimension $12\times 12$, are performed using 64 digits (octuple precision). The condition numbers for the mass and amplitude matrices (37) and (41) turn out to be equal to $\kappa(\overline{M})\approx 10^{31}$ and $\kappa(\overline{A})\sim 3\cdot 10^{5}$, respectively. Thus, with the octuple precision the numerical inversion of the mass and amplitude matrices is performed accurately.

We now apply the ODE algorithm [see Eqs. (36-42) of Section II.1] using the fake correlator $C_{n}^{(0)}$ and its derivatives $C_{n}^{(2k)}$ up to $k=N=12$, i.e. making the explicit use of the information on the number of exponential signals present in Eq. (68). The case in which, besides the correlator $C_{n}^{(0)}$, the derivatives $C_{n}^{(2k)}$ up to $k=N_{ODE}$ with $N_{ODE}\neq N$ are considered, is postponed to Sections III.2 and III.3. We remind that the definitions (37-38) and (41-42) correspond to the use of the values of the correlator $C_{n}^{(0)}$ in the full range $n=[1,\,N_{T}]$.

The ODE algorithm provides the values $X_{i}^{(ODE)}\equiv\{aM_{i}^{(ODE)},~{}A_{i}^{(ODE)}\}$ for the quantities $X_{i}\equiv\{aM_{i},~{}A_{i}\}$ and through the jackknife procedure the ODE errors $\delta X_{i}^{(ODE)}\equiv\{\delta(aM_{i}^{(ODE)}),~{}\delta A_{i}^{(ODE)}\}$ for the corresponding uncertainties $\delta X_{i}\equiv\{\delta(aM_{i}),~{}\delta A_{i}\}$ with $i=1,...,12$.

In the case at hand, for the relative deviations we get

[TABLE]

We have also considered the cases in which the uncertainties $\delta(aM_{i})$ and $\delta 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 $\Delta_{max}\sim 2\cdot 10^{-13}$ and $\delta\Delta_{max}\sim 3\cdot 10^{-11}$, while in the latter case we get $\Delta_{max}\sim 3\cdot 10^{-13}$ and $\delta\Delta_{max}\sim 2\cdot 10^{-9}$.

The above findings illustrate very clearly that the ODE algorithm is able not only to determine precisely the central values of masses and amplitudes, but also to detect accurately the fluctuations generated in the fake correlator by the uncertainties of masses and amplitudes. This ability is basically due to the huge value of the condition number $\kappa(\overline{M})$, i.e. to the closeness of the mass matrix to singularity. To deal with a huge value of $\kappa(\overline{M})$ is similar to the case of the study of a function $f(x)$ around a value $x=x_{0}$ where its derivative is huge. Small variations of $x$ around $x_{0}$ 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 1 are added in the fake correlator (68), both the condition numbers $\kappa(\overline{M})$ and $\kappa(\overline{A})$ as well as the maximum relative deviations $\Delta_{max}$ and $\delta\Delta_{max}$ increase quickly as shown in Table 2.

The quick rise of $\Delta_{max}$ and $\delta\Delta_{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_{n}^{(0)}$ 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 $\Delta_{max}$ and $\delta\Delta_{max}$). In Table 3 we have collected the values of $N_{max}$ found indicatively for different levels of the numerical precision of the correlator $C_{n}^{(0)}$ while keeping $\Delta_{max}<5\cdot 10^{-9}$ and $\delta\Delta_{max}<5\cdot 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 Eq. (68) we consider $N=4$ and a temporal extension equal to $N_{T}=T/a=48$. The values chosen for the amplitudes and the masses are collected in Table 4 together with their uncertainties $\delta A_{i}$ and $\delta(aM_{i})$. They have been chosen in such a way that at large time distances the correlator $C_{n}^{(0)}$ does exhibit a (short) plateau around half of the temporal extension of the lattice, as shown in the left panel of Fig. 2.

The plateau should in principle correspond to twice the amplitude $A_{1}$ of Table 4 (the factor of $2$ is due to the positive $t$-parity of the fake correlator), i.e. $2A_{1}=1.00(5)$. Instead, the average of the correlator in the plateau range $n=[22,\,26]$ provides the incorrect value $2A_{1}=1.16(5)$. Moreover, by adopting a simple fit based on two exponential signals in the larger range $n=[10,\,38]$ one gets an improved determination $2A_{1}=1.08(5)$, which however still differs from the input value. The temporal behavior of the effective mass $aM_{eff}^{(cosh)}$, given by Eq. (45) and shown in the right panel of 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 4.

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

[TABLE]

where $\mu_{i}$ is the multiplicity of the $i$-th exponential with $i=1,2,...\overline{N}$ and $\sum_{i=1}^{\overline{N}}\mu_{i}=N$. The values chosen for the masses $aM_{i}$, the amplitudes $A_{i\mu}$ and the multiplicities $\mu_{i}$ are collected in Table 5 together with the corresponding uncertainties $\delta(aM_{i})$ and $\delta A_{i\mu}$. 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.

The time dependencies of $C_{n}^{(2k)}$ for $k=0,...3$ and the one of the effective mass $aM_{eff}^{(cosh)}$ are shown in Fig. 3. It can be seen that at large time distances the effective mass (45) does not exhibit a plateau and, therefore, the lightest mass cannot be determined accurately by the standard effective mass procedure.

We now apply the ODE algorithm using all the derivatives $C_{n}^{(2k)}$ up to $k=N=12$. The condition numbers are found to be $\kappa(\overline{M})\approx 10^{28}$ and $\kappa(\overline{A})\sim 2\cdot 10^{12}$. Multiple roots are now present in the polynomial (58) due to the multiple poles in the fake correlator. The ODE algorithm determines accurately all the masses, amplitudes and multiplicities of Table 5. For the maximum relative deviations [see Eqs. (69-70)] we get $\Delta_{max}\sim 10^{-12}$ and $\delta\Delta_{max}\sim 2\cdot 10^{-10}$, which means that all the masses and amplitudes of Table 5 as well as the multiple pole structure of the correlator (71) are determined very precisely.

An interesting case is represented by the following fake correlator

[TABLE]

composed by six double poles with masses and amplitudes given in Table 6. Three signals ($m=1,3,5$) are non-oscillating double poles (i.e. the corresponding masses are real), while the other three signals ($m=2,4,6$) are oscillating double poles (i.e. the corresponding masses contain an imaginary part equal to $-\pi$). The numerical precision for generating the fake correlator (72) is again 32 digits (quadruple precision).

The application of the ODE algorithm (with 64 digits of precision) to the fake correlator $C_{n}^{(0)}$ and its derivatives $C_{n}^{(2k)}$ with $k=1,...,12$ is successful. All the six double poles are correctly found: three roots $\widetilde{z}_{m}\equiv 2\left[\mbox{cosh}(a\widetilde{M}_{m})-1\right]$ are positive (corresponding to real masses) and the other three roots are less than $-4$ (corresponding to masses with an imaginary part equal to $-\pi$). The condition numbers of the mass and amplitude matrices are $\sim 5\cdot 10^{30}$ and $\sim 3\cdot 10^{10}$, respectively. The accuracy of the ODE results for all the masses, amplitudes and their uncertainties is better than $\sim 1$ ppb.

Results with the same quality can be obtained in the case of multiple correlators (see Section II.4), 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) Blossier:2009kd . 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.

III.2 ODE analyses with $N_{ODE}<N$

In this Section we address the case in which the total number of derivatives $C_{n}^{(2k)}$ 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_{n}^{(0)}$. For the latter we use hereafter Eq. (68) with masses and amplitudes given in Table 1 for $N=12$.

We apply the ODE algorithm using $C_{n}^{(0)}$ and its derivatives $C_{n}^{(2k)}$ for $k=1,...,N_{ODE}$ with $N_{ODE}<N$, so that the mass and amplitude matrices have now dimension $N_{ODE}\times N_{ODE}$ . Furthermore, we consider that the range of analysis may not include all the values of the time distance at which the correlator $C_{n}^{(0)}$ is known, i.e. we may use a range $n=[n_{0}+1,\,N_{T}-n_{0}]$, where the integer $n_{0}$ may be larger than [math], at variance with what done in the previous Sections.

We write down explicitly the main ingredients of the ODE algorithm to take into account the use of a limited range of analysis $[n_{0}+1,\,N_{T}-n_{0}]$. The mass matrix (37) is replaced by

[TABLE]

with $k,k^{\prime}=0,1,...,(N_{ODE}-1)$ and the vector (38) is now given by

[TABLE]

The solution of the ODE conditions

[TABLE]

provides the coefficients $\overline{x}_{k}$, which are used to construct the polynomial of degree $N_{ODE}$

[TABLE]

having its roots at $z=\overline{z}_{j}^{ODE}\equiv 2\left[\mbox{cosh}(aM_{j}^{ODE})-1\right]$. Once the $N_{ODE}$ masses $aM_{j}^{ODE}$ have been determined, the corresponding amplitudes $A_{j}^{ODE}$ can be obtained by solving the linear system of equations

[TABLE]

where

[TABLE]

with333Eq. (80) holds for signals with real mass and multiplicity equal to $1$, as assumed in the fake data (68). Later on we generalize Eq. (80) to the case of imaginary ODE masses [see Eqs. (81-82)] and multiple poles [see Eqs. (84-85)].

[TABLE]

Generally speaking we expect that the ODE masses $aM_{j}^{ODE}$ represent an approximation of the lighter $N_{ODE}$ masses $aM_{j}$ included in the correlator $C_{n}^{(0)}$. The signals with masses $aM_{i}$ above $aM_{N_{ODE}}$ (i.e. for $i=[N_{ODE}+1,\,N]$) certainly affect the elements of both the mass matrix $\overline{M}$ and the vector $\overline{V}$. A kind of interaction among the roots may be generated, which moves the location of the ODE roots $\overline{z}_{j}^{ODE}$ of the polynomial (76) away from the exact ones $\overline{z}_{j}=2\left[\mbox{cosh}(aM_{j})-1\right]$, likely to higher values. The effect of the above interaction depends not only on the signals having the higher masses, but also on the range of the analysis, i.e. on the value of $n_{0}$. The difference between $aM_{j}^{ODE}$ and $aM_{j}$ is expected to decrease as $n_{0}$ increases, since the impact of the signals with higher masses decreases as the time distance increases.

This is indeed the phenomenology we observe. In Tables 7 and 8 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 there are masses with a non-vanishing imaginary part (namely $j=5$ and $j=6$). They are complex conjugate, since the coefficients of the polynomial (76) are real. Their presence is a signature that in the fake data of the correlator $C_{n}^{(0)}$ there are more than eight exponential signals. In order to calculate the amplitudes $A_{j}^{ODE}$ we have to take into account the presence of conjugate masses, which appear in pairs: one with a positive value of the imaginary part and the other with an opposite (negative) value of the imaginary part. Thus, we replace Eq. (80) with

[TABLE]

where444Eq. (81) holds for signals with multiplicity equal to $1$. Later on, Eq. (81) will be generalized to the case of multiple poles [see Eqs. (84-85)].

[TABLE]

In the limit ${\rm Im(M_{j}^{ODE})\to 0}$ the two exponential signals become a double pole and, indeed, Eq. (82) provides the appropriate limits $s_{n}^{j}\to 1$ for ${\rm Im(M_{j}^{ODE})\to 0^{+}}$ and $s_{n}^{j}\to n$ for ${\rm Im(M_{j}^{ODE})\to 0^{-}}$.

The lighter four exponential signals are reproduced within the uncertainties when the range of the analysis is equal to $n=[5,\,92]$, and quite precisely already when $n=[10,\,87]$. On the contrary the heavier four signals are not reproduced at all for $n=[5,\,92]$ and only the mass of the fifth signal matches the exact one within the uncertainty for $n=[10,\,87]$.

Though not all the ODE masses and amplitudes reproduce the input values, the fake data of the correlator $C_{n}^{(0)}$ are reasonably reproduced by its ODE representation

[TABLE]

We find that the absolute value of the residue $[C_{n}^{(0)}-C_{n}^{(ODE)}]$ does not exceed $2\cdot 10^{-3}\,\sigma_{n}^{(0)}$ and $10^{-5}\,\sigma_{n}^{(0)}$ for the analyses in the ranges $n=[5,\,92]$ and $n=[10,\,87]$, respectively. This means that small values of the residue do not guarantee that masses and amplitudes are properly reproduced.

The remarkable stability of the masses and amplitudes of the lighter four states, determined by the ODE algorithm with $N_{ODE}=8$, for various ranges of the analysis is also illustrated in Fig. 4, where the results for the ratios $M_{j}^{ODE}/M_{j}$ and $A_{j}^{ODE}/A_{j}$ are shown.

On the contrary, an approximate stability is found for the fifth and sixth states only at large values of $n_{0}$, as shown in Fig. 5. A very bad convergence of the masses and amplitudes of the heavier two ODE states ($j=7$ and $j=8$) is clearly visible.

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_{n}^{(2k)}$ is more sensitive to the signals with higher masses as the order $k$ increases (cf., e.g., the factor $(\overline{z}_{i})^{k}$ in Eq. (34)). This feature is illustrated in Fig. 6, where for a given choice of the analysis range ($n=[20,\,76]$) the convergence for the masses and amplitudes of the ODE states with $j=5,..,8$ toward the input values is clearly visible.

For sake of completeness we now generalize the structure of the ODE representation (83) to the case of multiple exponential signals having masses $M_{j}^{ODE}$ with non-vanishing imaginary parts and multiplicity $\mu_{j}$ greater than $1$. One has

[TABLE]

where $\sum_{j=1}^{\overline{N}_{ODE}}\mu_{j}=N_{ODE}$ and

[TABLE]

with $s_{n}^{j}$ given in Eq. (82) and $(-)^{p}$ being the $t$-parity of the correlator.

Finally, in the ODE procedure we introduce two tolerance factors $\delta_{R}$ and $\delta_{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

[TABLE]

the two signals are replaced by a signal having as real part $[\mu_{j}a{\rm Re}(M_{j}^{ODE})+\mu_{j^{\prime}}a{\rm Re}(M_{j^{\prime}}^{ODE})]/(\mu_{j}+\mu_{j^{\prime}})$ and multiplicity $\mu_{j}+\mu_{j^{\prime}}$. When the imaginary part of the mass of a signal is small enough, namely

[TABLE]

it is put to zero. By using the above tolerance factors the relative error in the ODE representation (84)) is of the order ${\cal{O}}(\delta^{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.

III.3 ODE analyses with $N_{ODE}>N$

In Section III.1 we have illustrated the nice performance of the ODE algorithm applied to the case in which the number of derivatives of the correlator included in the construction of the mass matrix is equal to the number of exponential signals included in the fake correlator, i.e. $N_{ODE}=N$. 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 3).

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 $\overline{P}_{N_{ODE}}(z)$ (see Eq. (76)) is larger than the number $N$ of signals present in the correlator $C_{n}^{(0)}$.

Let us consider again the fake data for the correlator (68) corresponding only to the first eight masses and amplitudes given in Table 1 with $T/a=N_{T}=96$. When $N_{ODE}=N=8$ the ODE algorithm determines precisely the input masses and amplitudes (together with their statistical fluctuations). The absolute residue $\left|C_{n}^{(0)}-C_{n}^{(ODE)}\right|$ between the input correlator and its ODE representation is numerically very small ($<10^{-20}$), and it corresponds to the level of the numerical rounding of the whole procedure.

We now apply the ODE algorithm using $N_{ODE}=10,12,14,16$. The condition number of the mass matrix quickly increases up to $\approx 10^{70}$ for $N_{ODE}=16$. Therefore we adopt 96 digits for the internal ODE computations. For all values $N_{ODE}>8$ the eight roots of the fake correlator are properly determined, but also additional ($N_{ODE}-8$) roots are found. The interesting feature is that the values of $\overline{z}_{j}^{ODE}\equiv 2\left[\mbox{cosh}(aM_{j}^{ODE})-1\right]$ for $j>8$ are real and in the range $(-4,\,0)$, which in turn mean imaginary values of the corresponding ODE masses $aM_{j}^{ODE}$, as shown in Table 9.

The extra ($N_{ODE}-8$) roots correspond to noisy signals, that the ODE algorithm includes for trying to take into account the small residual terms. Note that the extra masses of Table 9 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_{j}^{ODE}$ 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 3 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_{PS}(t)$ defined as

[TABLE]

where $P_{5}(x)=\overline{q}_{2}(x)\gamma_{5}q_{1}(x)$ is a local interpolating field that creates in $x$ a pseudoscalar (PS) mesons made of two valence quarks $\bar{q}_{1}$ and $q_{2}$ with masses $m_{1}$ and $m_{2}$. At large time distances one has

[TABLE]

so that the meson mass $M_{PS}$ and the matrix element $\mathcal{Z}_{PS}=\langle PS|\overline{q}_{2}\gamma_{5}q_{1}|0\rangle$ 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 McNeile:2006bz , 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 Ref. Carrasco:2014cwa ), we have selected the case of the gauge ensemble B25.32 corresponding to a lattice volume $V\times T=32^{3}\times 64\,a^{4}$, to a lattice spacing $a\simeq 0.082$ fm with a light-quark mass equal to $m_{\ell}\simeq 12$ MeV (in the $\overline{\rm MS}(2~{}\mbox{GeV})$ scheme) and a strange and charm masses close to their physical values. The number of independent gauge configurations employed is 150. The correlator (88) has been calculated using 160 stochastic sources (diagonal in the spin variable and dense in the color one) per gauge configuration for valence quark masses equal to $m_{1}=m_{2}=m_{\ell}$, corresponding to a simulated pion mass of $\simeq 260$ MeV (see Ref. Giusti:2018mdh ). The relative statistical error $\sigma_{PS}(t)/C_{PS}(t)$ does not exceed $\simeq 1.4\%$.

The ETMC correlator (88) and its second and fourth derivatives are shown in the left panel of Fig. 7. With respect to the case of fake data (see the left panel of Fig. 1) the fourth derivative (as well as higher ones) is manifestly noisy at large time distances.

We have applied the ODE algorithm using values of $N_{ODE}$ from $4$ up to $10$ in the full range of values of $n$, i.e. $n=[1,\,64]$. In the case $N_{ODE}=4$ we find three exponential signals, which are non-noisy and non-oscillating signals (i.e. $\overline{z}_{j}^{ODE}>0$), and $1$ noisy signal of the type described in the previous Section (i.e. $-4<\overline{z}_{j}^{ODE}<0$).

In what follows we refer to the non-noisy states as the physical states.

When $N_{ODE}>4$ we find always four physical signals and $N_{ODE}-4$ noisy states. For all values of $N_{ODE}$ physical and noisy states have multiplicity equal to $1$. Thus, for the correlator (88) the ODE representation corresponding only to the physical states is given by

[TABLE]

where $N_{phys}$ indicates the number of physical states. The difference $C_{PS}(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$. The noisy states interfere clearly with the detection of excited states and the number of physical exponential signals, that can be extracted by the ODE method, is therefore limited by the level of the noise. We speculate that the observed noise can be reduced by increasing the number of gauge configurations used. This point requires, however, further dedicated investigations.

The results for the physical signals obtained with $N_{ODE}=4$, $5$, $6$, $8$ and $10$ are collected in Table 10.

A good convergence for the ground state and progressively for the excited states is observed as $N_{ODE}$ increases. For $N_{ODE}\geq 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}\geq 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

[TABLE]

At enough large time distances both definitions provide the mass of the ground-state. However, at lower values of $n$ they differ due the different impact of the excited states.

A single exponential fit in a chosen range $n=[n_{min},\,n_{max}]$ provides values for $aM_{1}$ and $A_{1}$, which in general may depend on the choice of $n_{min}$ and $n_{max}$. Such a dependence may represent a possible source of a systematic uncertainty. In the case at hand a single exponential fit in the range $[18,\,30]$ yields $aM_{1}=0.10699~{}(37)$ and $A_{1}=0.08025~{}(55)$ in nice agreement with the corresponding ODE results of Table 10. The dependence of the chosen range turns out to be a sub-leading effect with respect to the statistical error.

We stress that the ODE procedure allows to extract the ground-state signal taking into account the presence of excited states without the need of finding a range of values of $n$, where the ground-state dominates. The systematic uncertainty related to the choice of such a range (i.e. to the contamination of excited states) is therefore avoided adopting the ODE procedure. Moreover, important pieces of information on the excited states can be obtained for spectroscopic studies.

We now move to the heavier sector of charmonium. We select the gauge ensemble D20.48 Carrasco:2014cwa corresponding to a lattice volume $V\times T=48^{3}\times 96a^{4}$ and to a lattice spacing $a\simeq 0.062$ fm, and choose the valence quark masses equal to $m_{1}=m_{2}\simeq 1.18$ GeV (in the $\overline{\rm MS}(2~{}\mbox{GeV})$ scheme), i.e. close to the physical charm quark mass Carrasco:2014cwa . The number of independent gauge configurations employed is 100. We consider both the PS correlator (88) and the vector (V) one given by

[TABLE]

where $V_{i}(x)=\overline{q}_{2}(x)\gamma_{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 $\sim 0.4\%$ and $\sim 0.7\%$ for the PS and V correlators, respectively.

We apply the ODE algorithm with $N_{ODE}\geq 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 11 and 12 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. In the PS case a single exponential fit in the range $n=[26,\,44]$ provides the values $aM_{1}=0.95175~{}(14)$ and $A_{1}=0.06426~{}(25)$, which agree well with the corresponding ODE results shown in Table 11.

In the V case the plateau of the effective mass is of lower quality and the dependence of the single exponential fit on the chosen range $[n_{min},\,n_{max}]$ is more relevant, as it is shown in Table 13.

The corresponding systematic uncertainty turns out to be not negligible with respect to the statistical error. Note also that the results of the single exponential fit lie below the corresponding ODE results, shown in Table 12, by approximately $2$ standard deviations.

V Use of the ODE algorithm together with other techniques

The results presented in the previous Sections show that the use of the ODE algorithm allows to extract multiple exponential signals from the temporal dependence of correlation functions without the need of any further input other than the correlator itself. In particular, the ODE method is able to detect the multiplicities of the signals as well as the presence of possible oscillating signals. This represents a relevant piece of information about the specific structure of the temporal dependence of the correlator.

In this way even a quite involved representation, like the one given by Eq. (84), can be explicitly written down and the ODE algorithm provides values for both masses and amplitudes as well as their uncertainties. All that can be used as a starting point for the the application of other techniques suitable for refining the fitting procedure of the temporal dependence of the correlator.

The simplest choice is to adopt a nonlinear least-squares minimizer starting from the ODE solution corresponding to the physical states. In other words, the ODE physical states provide the specific structure of the temporal dependence of the correlator, where masses and amplitudes can be used as free parameters to be varied starting from the values obtained by the ODE algorithm and to be determined by minimizing a $\chi^{2}$-variable. If the quality of the ODE representation of the correlator is adequate, then the subsequent application of the least-squares minimizer will either confirm the ODE values within the uncertainties or possibly refine the ODE solution.

Note that the ODE method does not correspond to the minimization of a unique $\chi^{2}$-variable. Indeed, as pointed out in Section II, the inversion of the mass matrix is equivalent to minimize the variable $\chi_{M}^{2}$ defined by Eq. (21), while, once the masses are given, the determination of the amplitudes corresponds to the minimization of the $\chi^{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 $\chi^{2}$-minimization procedure is expected to be sensitive also to the amplitudes of the noise. In this way the combined ODE plus $\chi^{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 $\chi^{2}$-minimization procedure. An interesting feature is that no priors on the masses and amplitudes are required in the $\chi^{2}$ minimization in order to obtain stable results. In the constrained curve fitting method of Ref. Lepage:2001ym 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. Giusti:2018mdh . 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. mpfun . Moreover, the root finding has been carried out accurately using the open-source software package MPSolve MPSolve .

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 detected by the ODE algorithm when analyzing lattice correlators, requires dedicated numerical investigations, that we are planning. At the same time applications of the ODE algorithm to the analysis of 2-point and 3-point correlation functions evaluated on the lattice are in progress.

Acknowledgments

We gratefully acknowledge V. Lubicz and F. Sanfilippo for fruitful discussions and D. Giusti for a careful reading of the manuscript. S.S. warmly thanks D.H. Bailey for his valuable assistance in the use of the MPFUN2015 software package.

Appendix A Filtering of signals

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

[TABLE]

The first step of the ODE algorithm is to calculate iteratively the derivatives

[TABLE]

for $k=1,...,N$ and construct the mass matrix (37) and the vector (38). Then, the linear system of equations (36) can be solved numerically to obtain the coefficients $x_{k}$, which define the polynomial

[TABLE]

having $N$ roots located at $z_{i}=2\left[\mbox{cosh}(aM_{i})-1\right]$.

Now we want to subtract a subset of exponential signals, for instance $M$ exponentials having the masses $M_{i}$ with $i=1,...,M$. First of all, we rewrite the polynomial (95) as

[TABLE]

where

[TABLE]

Since the roots of $Q_{M}(z)$ are known, it is easy to calculate explicitly the coefficients $y_{k}$ of the polynomial $Q_{M}(z)$, viz.

[TABLE]

with

[TABLE]

We now construct the modified correlator

[TABLE]

which implies

[TABLE]

Since by definition $Q_{M}(z_{i})=0$ for $i=1,...,M$, the modified correlator does not contain any more the unwanted exponential signals, i.e. those corresponding to the roots $z_{i}$ with $i=1,...,M$, while it contains only the remaining $N-M$ exponential signals

[TABLE]

with

[TABLE]

We now apply the ODE algorithm to the modified correlator $\widetilde{C}_{n}^{(0)}$, check that it has $N-M$ roots at $z_{i}$ with $i=M+1,...,N$ and obtain the amplitudes $\widetilde{A}_{i}$. Finally we reconstruct the original amplitudes $A_{i}$ using Eq. (104).

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 (93). Then, the filtering of the unwanted states can proceed by constructing the modified correlator (101) and by applying to it the ODE algorithm.