If one is interested in more than the ground state in a channel with given quantum numbers, the method of choice is the so-called generalised eigenvalue method (GEVM) based on a generalised eigenvalue problem (GEVP). The GEVM can also be used to extract ground state properties if there are large excited state contaminations such that from a single correlation function estimates are unreliable. The references to read are (Michael and Teasdale 1983; Lüscher and Wolff 1990; Blossier et al. 2009 ; Fischer et al. 2020)
Consider correlator matrices in the spectral representation of the form \[\begin{equation} C_{ij}(t)\ =\ \langle \hat O_i^{}(t')\ \hat O_j^\dagger(t'+t)\rangle\ =\ \sum_{k=0}^\infty e^{-E_k t} \psi_{ki}^*\psi_{kj}\,, \end{equation}\] with energy levels \(E_k > 0\) and \(E_{k+1} > E_k\) for all values of \(k\). The \(\psi_{ki}=\langle 0|\hat O_i|k\rangle\) are matrix elements of \(n\) suitably chosen (local) operators \(\hat O_i\) with \(i=0, ..., n-1\). Then, the eigenvalues or so-called principal correlators \(\lambda^0(t, t_0)\) of the generalised eigenvalue problem (GEVP) \[\begin{equation} C(t)\, v_k(t, t_0)\ =\ \lambda^0_k(t, t_0)\, C(t_0)\, v_k(t, t_0)\,, \end{equation}\] can be shown to read \[\begin{equation} \label{eq:lambdaGEVP} \lambda^0_k(t, t_0)\ =\ e^{-E_k(t - t_0)} \end{equation}\] for \(t_0\) fixed and \(t\to\infty\). Clearly, the correlator matrix \(C(t)\) will for every practical application always be square but finite with dimension \(n\). This will induce corrections to the form of \(\lambda\). The corresponding corrections were derived in Ref. (Lüscher and Wolff 1990; Blossier et al. 2009) and read to leading order \[\begin{equation} \label{eq:t0corr} \lambda_k(t, t_0) = b_k \lambda^0_k(1+\mathcal{O}(e^{-\Delta E_k t})) \end{equation}\] with \(b_k>0\) and \[\begin{equation} \Delta E_k\ =\ \min_{l\neq k}|E_l -E_k|\,. \end{equation}\] Most notably, the principal correlators \(\lambda_k(t_0, t)\) are at fixed \(t_0\) again a sum of exponentials. As was shown in Ref. Blossier et al. (2009), for \(t_0>t/2\) the leading corrections are different, namely of order \[\begin{equation} \label{eq:t0corr2} \exp[-(E_n - E_k)t]\,. \end{equation}\]
In practice, we need to solve the GEVP \[\begin{equation} C(t)\, v_k(t, t_0)\ =\ \lambda^0_k(t, t_0)\, C(t_0)\, v_k(t, t_0)\,, \end{equation}\] for each \(t\) and \(t_0\). We will assume that the correlator matrix is symmetric positive definite. You should always symmetrise your matrix, even if you know that it is symmetric in theory, to avoid numerical problems.
Now, compute the Cholesky decomposition of \(C(t_0)=L^t\cdot L\) and re-write \[ (L^{-1})^t \cdot C(t) \cdot L^{-1}\, L\,v_k\ =\ \lambda^0_k L\, v_k\,. \] Wet \(w_k = L v_k\) to obtain an ordinary eigenvalue problem \[ M(t, t_0)\, w_k = \lambda^0_k w_k\,. \] with \(M(t, t_0)=(L^{-1})^t \cdot C(t) \cdot L^{-1}\) preserving the symmetry.
After this step you will have \(\lambda^0(t, t_0)\) for certain combinations of \(t\)- and \(t_0\)-values. The next step consists of sorting the eigenvalues (and the corresponding eigenvectors) to obtain the consistently for each \(k\). The easiest way to do this is sorting by magnitude. Note that this might lead to incorrect matching in the case of noise. For more advanced ways to match see Ref. Fischer et al. (2020).
The goal of this exercise is to implement the GEVP for an artificial dataset without noise. Thus, no error estimate will be neccessary.
First implement the GEVP described above for a fixed pair of matrices \(C(t)\) and \(C(t_0)\). Include a step that makes sure matrices are positive definite. Note that it is not needed to assume an order between \(t\) and \(t_0\). Also include the case \(t=t_0\).
with \(t\) running. For the analysis it will be useful to define appropriate effective masses.
Analyse the artificial data and determine the corresponding energy levels. Verify the formulae from above for the corrections due to residual excited states.
Blossier, Benoit, Michele Della Morte, Georg von Hippel, Tereza Mendes, and Rainer Sommer. 2009. “On the generalized eigenvalue method for energies and matrix elements in lattice field theory.” JHEP 04: 094. https://doi.org/10.1088/1126-6708/2009/04/094.
Fischer, Matthias, Bartosz Kostrzewa, Johann Ostmeyer, Konstantin Ottnad, Martin Ueding, and Carsten Urbach. 2020. “On the generalised eigenvalue method and its relation to Prony and generalised pencil of function methods.” Eur. Phys. J. A 56 (8): 206. https://doi.org/10.1140/epja/s10050-020-00205-w.
Lüscher, Martin, and Ulli Wolff. 1990. “How to Calculate the Elastic Scattering Matrix in Two-dimensional Quantum Field Theories by Numerical Simulation.” Nucl. Phys. B339: 222–52. https://doi.org/10.1016/0550-3213(90)90540-T.
Michael, Christopher, and I. Teasdale. 1983. “Extracting Glueball Masses from Lattice QCD.” Nucl. Phys. B215: 433–46. https://doi.org/10.1016/0550-3213(83)90674-0.