\chapter{Exploratory techniques}
\label{cha:Exploratory}
 Most of the works contributing to the muon $g-2$ puzzle only give results for certain windows. It was shown that a large amount of the uncertainties are related to the bad signal-to-noise ratio of the vector correlator at large time distances. While in the fist high-precision estimation of $a_\mu^\mathrm{LO-HVP}$~\cite{Borsanyi_2021} the bounding method was used (which uses estimates for upper and lower bound the vector propagator in the critical regime), the idea came up to reconstruct the long-distance behavior by a spectral decomposition of the correlation function at shorter distances~\cite{Lahert:2021xxu,Lahert:2023ore,Djukanovic:2024cmq,RBC:2024fic}. Nevertheless, the by far best reduction of statistical uncertainties in the long-distance regime could be made by including the data-driven approach for only a small part the long-distance window while using the lattice simulations for the bulk of $a_\mu^\mathrm{LO-HVP}$. This approach was first published by the BMW collaboration in 2024~\cite{Boccaletti:2024guq}. Before it could shown that the discrepancies of data-driven and lattice approach arise from the shorter distances~\cite{Davier:2023cyp}. Even though we have a very precise methods for the tail distributions by now, it is still interesting to reconstruct the vector propagator for large distances to have an independent check of the data-driven approach.\\
\\
Since there is a strong decay of the vector meson ($\rho$) into two pions, the vector correlator has contributions of the form $\exp(-E_{2\pi} t)$\footnote{$E_{2\pi}$ denotes the energy of a two-pion system, thus, $E_{2\pi} \geq 2 M_\pi$.}. If there the are two-pion states with an energy below the mass of the $\rho$-meson (which is actually the case), the long-distance regime is dominated by the two-pion states. So we plan to extract the two-pion contributions from the vector correlator using a GEVP (see Section~\ref{sec:GEVP}) that includes the vector meson and two-pion states at short distances. The extracted parameters can then be used to reconstruct the correlator at long distances.\\
This chapter is divided into the following sections. First, I present my mathematical consideration to construct two-pion states with the same quantum numbers as the vector meson using staggered fermions. Next, I will present the ideas of the implementation. In particular I will construct a method for an efficient computation of connected two-pion diagrams. The content of the following section has been presented on the LATICE2023 and LATTICE2024 conferences and published in the corresponding Proceedings of Science (PoS)~\cite{Frech:2023bzy,Frech:2025fjp}.\\



\section{Two pion contribution to vector mesons}
\label{sec:pipi_math}
The staggered symmetry group on a fixed time slice is given by
\begin{align}
  \mathbb{Z}_N^3 \rtimes \left[ \Gamma^{4,1} \rtimes W_3\right],
\end{align}
with $\rtimes $ denoting the semi-direct product. The first factor describes the translations on the three dimensional torus, the second one the charge conjugation and taste transformations and the last factor the rotations and inversions of three dimensional cube. More details on the symmetry group of the staggered action is given in Section~\ref{sec:stagg_symm}. In this study a vector meson at rest, trivial in taste space, is considered. Therefore, the representation of the vector state is given by
\begin{align}
C \otimes V,
\end{align}
where $C$ is the non-trivial ($\mathbb{Z}_2$) representation of the charge conjugation\footnote{Later we will include $\pi\pi$ states with negative charge conjugacy ($I = 1$ combination). Therefore, we can ignore the charge conjugation issues in the mathematical derivation.} and $V$ is the vector representation of the Würfel group, i.e., the defining representation of the orthogonal group $O(3)$ restricted on $W_3$.\\
The pion states are somewhat more involved since they are not trivial in position and taste. Indeed, the representation of resting Goldstone-pions is trivial. However, when considering two-pion states taste momentum structure will be more involved. Basis vectors of the pion representation are given by
\begin{align}
\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm},
\end{align}
where the 4th component of the taste is separated from the others, because this component is not affected by the spatial rotations. These states transform under the generators of the staggered symmetry group in the following way
\begingroup
\allowdisplaybreaks
\begin{align}
S_{i}\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm} &= e^{-ip_i}\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm}\\
\Xi_{\mu}\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm} &= (-)^{\xi_{\mu}}\ket{\vec p,\vec \xi, \xi_{4}}_{\pi}\\
R\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm} &= \det(R)\ket{R\vec p,R\vec \xi, \xi_{4}}_{\pi^\pm}\\
C\ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\pm} &= \ket{\vec p,\vec \xi, \xi_{4}}_{\pi^\mp}.
\end{align}
\endgroup
This structure can be derived using the little group procedure mentioned in the Appendix~\ref{sec:wigner}.\\
A two-pion state with $I = 1$ is given (up to normalization) by the wedge product $a \wedge b = \frac{1}{\sqrt{2}}(a\otimes b - b\otimes a)$ of two single-pion states with charge $+$ resp. $-$. So the group representation is generated by
\begin{align}
S_{i}\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right) &= e^{-i(p^a_i + p^b_i)}\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right)\\
\Xi_{\mu}\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right) &= (-)^{\xi^a_{\mu} + \xi^{b}_{\mu}}\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right)\\
R\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right) &= \left(\ket{R\vec p^a,R\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{R\vec p^b,R\vec \xi^b, \xi_{4}^b}_{\pi^-}\right)\\
C\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right) &= -\left(\ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-}\right).
\end{align}
Here, it is used that $\det(R)^2 = 1$ for all elements of $W_3$. For simplicity we will denote these states as
\begin{align}
 \ket{\vec p^a,\vec \xi^a, \xi_{4}^a}_{\pi^+}\wedge \ket{\vec p^b,\vec \xi^b, \xi_{4}^b}_{\pi^-} \eqqcolon \ket{\vec p^a,\vec \xi^a, \xi_{4}^a;\vec p^b,\vec \xi^b, \xi_{4}^b}.
\end{align}
Of course, not all two-pion states have the same overlap to all vector meson states. Quantum numbers have to be conserved as well as some of the states appear with a different multiplicity. Mathematically spoken, it behaves the following way: The vector meson representation of the Würfel group is irreducible while the two-pion representations are reducible. Therefore, the two-pion states can be decomposed into irreducible representations, from which we expect some to be the vector meson irrep. So we can construct vector-meson like states from the two-pion states with the same quantum numbers according to
\begin{align}
\label{eq:meson_state_decomp}
\begin{aligned}
  & \ket{\vec q = 0, \zeta_\mu = 0; s,m,\alpha}= 
   \sum_{(\vec p^a,\xi_{\mu}^a,\vec p^b,\xi_{\mu}^b)\in s} &\ket{\vec p^a,\vec \xi^a, \xi_{4}^a;\vec p^b,\vec \xi^b, \xi_{4}^b} \braket{\vec p^a,\vec \xi^a, \xi_{4}^a;\vec p^b,\vec \xi^b, \xi_{4}^b}{\vec q = 0, \zeta_\mu = 0; s,m,\alpha},
\end{aligned}
  \end{align}
  where $\vec q$ denotes the momentum and $\zeta_{\mu}$ the taste of the vector meson, $s$ labels the orbit $O_s$ of the two-pion states under $W_3$ rotations and spatial inversions, $m \in (1,...,a_s)$ labels the multiplicities with which the corresponding orbit occurs, and $\alpha = x,y,z$ represents the orientation of the reconstructed vector meson.
One should note that the vector-like states have a rotational symmetry and therefore, the final states only depend on the absolute value of the momentum and taste vectors (this is why they are classified into orbits) and on the $4^{th}$ component of the taste. Of course, the rotations affects the direction of the vector state but this is considered in further calculations.\\ 
The derivation is similar to Wigner's little group procedure mentioned in Section~\ref{sec:wigner} and in~\cite{Weinberg}. \\
So one has to figure out the multiplicities $a_s$ as well as the Clebsch-Gordan coefficients which are given defined as the scalar product from Equation~\eqref{eq:meson_state_decomp}. Computing the multiplicities is straightforward using basic essences of finite group theory. The character of a representation is given by the trace of a corresponding group element, which means in the case of a two-pion state
\begin{align}
\begin{aligned}
  &\chi^{2\pi}(\vec S,R,\Xi_{\mu}) = \\
\sum_{\vec p^a, \vec p^b, \xi_\mu^a, \xi_\mu^b} &\braket{\vec p^a,\vec \xi^a, \xi_{4}^a;\vec p^b,\vec \xi^b, \xi_{4}^b}{R\vec p^a,\vec R\xi^a, \xi_{4}^a;R\vec p^b,R\vec \xi^b, \xi_{4}^b}\,  e^{-i \vec S\cdot R(\vec p^a + \vec p^b) - i\pi \vec \Xi\cdot R(\vec \xi^a + \vec \xi^b) -i\pi\Xi_4(\xi^a_4 + \xi^b_4)}
\end{aligned}\\
\end{align}
and in the case of the vector meson
\begin{align}
\begin{aligned}
  \chi^{\rho}(R) = \Tr(R).
\end{aligned}
\end{align}
Then, the multiplicity of a given orbit is given by
\begin{align}
a_s &= \frac{1}{2^4\cdot 48\cdot L^3}\sum_{nR\Xi}\chi^\rho(R)^* \times \chi^{2\pi}(n,R,\Xi)\\
 &=\frac{1}{48}\sum_{R}\chi_\rho(R)^*\times \sum_{\vec p^b}\sum_{\vec p^a =R\vec p^a}\delta(\vec p^a+\vec p^b)\times \sum_{\vec\xi^b}\sum_{\vec \xi^a = R\vec\xi^a}\delta(\xi^a_{1,2,3,4} - \xi^b_{1,2,3,4}),
\end{align}
which can be computed very easily using a small script. Please note that the resulting formula is independent of the lattice volume as long as the cubic symmetry is valid.\\
When computing this coefficients one will notice that the only two-pion states that have overlap to the resting taste-singlet vector meson are those that have total momentum $0$ and constitute of two pions that have the same taste. Therefore, the multiplicities (and the CG-coefficients) can be labeled by a single momentum (which we call $\vec p$) and a single taste (which we call $\xi$)\footnote{Of course, $\vec p$ can be replaced by $-\vec p$. This results in an interchange of the two pions and produces only a global sign which is not relevant in the generalized eigenvalue problem.}. Independent of the choice of $\xi_4$ we get the coefficients shown in table \ref{tab:multies}. The multiplicities as well as the Clebsch-Gordan coefficient only depend on the orbit of the pion state under the action of the Würfel group. The orbits can be labeled by the absolute value of taste and momentum and from this point we will denote the two-pion states that have total momentum zero and that are taste-singlets in the bra-ket notation as
\begin{align} 
 \ket{\vec p,\vec \xi, \xi_{4};-\vec p,\vec \xi, \xi_{4}} \eqqcolon \ket{\vec p,\vec \xi,\xi_4; I = 1}.
\end{align}
Next one needs to compute the Clebsch-Gordan coefficients, that are given by the scalar product in Equation~\eqref{eq:meson_state_decomp}:
\begin{align}
  C^{s,m, \alpha}_{\vec p, \vec \xi, \xi_4; R} = \braket{R\vec p, R\vec \xi,\xi_4; I = 1}{\vec q = \vec 0, \zeta_\mu = 0; s, m,\alpha}.
\end{align}
The Clebsch-Gordon coefficients do not depend on $\xi_4$ and can of course be given in terms of momenta and tastes from the corresponding rotational orbit: 
\begin{align}
  C^{s,m, \alpha}_{\vec p, \vec \xi} = \braket{\vec p, \vec \xi,\xi_4, I = 1}{\vec q = \vec 0, \zeta_\mu = 0;s,m,\alpha}.
\end{align}
Since the symmetry group that is considered here is a finite group, we can use brute-force techniques~\cite{CG-Finite-Group} in order to compute the CG-coefficients. The application of this method for the two-pion-vector-meson overlap is explained very briefly in the following paragraph:\\
Let $T(R)$ be the vector representation of a Würfel group element $R$. In this case, all of those 48 matrices are real, orthogonal and three dimensional. For a given taste and momentum orbit $O$ and Rotation $R\in W_3$ we define $D^O(R)$ as a matrix of dimension $\vert O \vert$ fulfilling the following condition:
\begin{align}
\label{eq:twopionrep}
D^O(R)_{\vec p^a,\xi^a;\vec p^b,\xi^b} = \delta_{\vec p^a,R\vec p^b}\delta_{\xi^a,R\xi^b}\quad\forall (\vec p^a,\xi^a),\, (\vec p^b,\xi^b) \in O
\end{align}
One has to note that for the taste one does not need to distinguish between positive and negative number, i.e., the rotations $R$ act on the taste vector as elements of $O(\mathbb{Z}_2^3)$. For example a rotation be $\pi$ around the $z$-axis leaves the unit vector in $z$-unit vector invariant:
\begin{align}
\begin{pmatrix}
1&0&0\\0&-1&0\\0&0&-1
\end{pmatrix}
\begin{pmatrix}
0\\0\\1
\end{pmatrix}
\equiv
\begin{pmatrix}
0\\0\\1 
\end{pmatrix}\mod 2.
\end{align}
One might ask why the representation in Equation \eqref{eq:twopionrep} has no dependence on the shift or taste operation. This is because in this special case of two pion states the phases (signs) for shift (taste) operation cancel exactly. E.g. while the first pion gives a phase of $e^{i\vec p^a \cdot\vec s}$, the second one gives a phase of $e^{i\vec p^b \cdot\vec s} = e^{-i\vec p^a \cdot\vec s}$. Therefore, the whole computation becomes independent of the lattice size which reduces the computational effort by far. The whole computation can be done within a few seconds on a usual laptop PC.\\
We can compute the coefficients by normalizing the columns of
\begin{align}
\label{eq:CG}
\sum_R D^O(R)\cdot A \cdot T(R)\, ,
\end{align}
where $A$ is a matrix with general coefficients and the size of $A$ is chosen such that the multiplications are well-defined.
Most of the general entries of $A$ will cancel after this procedure but the remaining ones will span a vector space whose dimension is the multiplicity we computed before\footnote{In principle we do not need to compute the multiplicities before, but it is very straightforward and some kind of crosscheck for the computation of the coefficients.}. Choosing a certain basis here will give the coefficients for the different copies of the vector space. The final matrix will have three columns each representing one component of the vector current. The entries in each of the rows will denote the Clebsch-Gordan coefficient of the corresponding state. For example let us consider the $\vert \vec p \vert = 1$, $\vert \vec \xi \vert = 3$ and $\xi_4 = 1$ orbit with the ordered basis
\begin{align}
\ket{\vec p = \vec e_x, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1}, \ket{\vec p = -\vec e_x, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1},\\
\ket{\vec p = \vec e_y, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1}, \ket{\vec p = -\vec e_y, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1},\\
\ket{\vec p = \vec e_z, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1}, \ket{\vec p = -\vec e_z, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1}.
\end{align}
The resulting matrix (after normalizing the columns) of Equation~\eqref{eq:CG} will be
\begin{align}
\frac{1}{\sqrt{2}}
\begin{pmatrix}
1&0&0\\
-1&0&0\\
0&1&0\\
0&-1&0\\
0&0&1\\
0&0&-1
\end{pmatrix},
\end{align}
which means that the $i$th component of the vector overlap with this orbit is given by
\begin{align}
\frac{1}{2}\left(\ket{\vec p = \vec e_i, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1} -\ket{\vec p = -\vec e_i, \vec \xi = \vec e_x + \vec e_y + \vec e_z, \xi_4 = 1}\right).
\end{align}

For states with a multiplicity larger than one (e.g. the same pion states with $\vert \vec \xi \vert = 2$), one can construct more independent linear combinations of the overlapping states.

Next we will present the CG coefficients for the taste and momentum orbits up to $\vec p^2 = 4$. Due to symmetries, there is no distinction between $\vec p^2 = 1$ and $\vec p^2 = 4$. The unit vectors in momentum space are denoted with $\vec e_i$, while $f_i$ represents unit vectors in taste space, whenever $\vec p^2 = 1$ or $\vec \xi^2 = 1$. If the orbits are of the type $\vec \xi^2 = 2$, $f_i$ denotes a vector with a zero in the $i$th component and a one in the other ones. $\vec 0$ denotes the zero vector in the case $\vec \xi^2 = 0$ and $(1,1,1)^T$ when $\vec \xi^2 = 3$. Here, $\lambda$, $\mu$, and $\nu$ refer to $1$ and $-1$, respectively. The Clebsch-Gordan coefficients are shown in Table \ref{tab:CG}.
\begin{table}[]
  \centering
  \begin{tabular}{c||c|c|c|c|c}
              &$\vert \vec p \vert^2 = 0$&$\vert \vec p \vert^2 = 1$&$\vert \vec p \vert^2 = 2$&$\vert \vec p \vert^2 = 3$&$\vert \vec p \vert^2 = 4$\\
              \hline\hline
     $\vert \vec \xi\vert^2 = 0$  & 0&1&1&1&1\\
       $\vert \vec \xi\vert^2 = 1$& 0&2&3&2&2\\
       $\vert \vec \xi\vert^2 = 2$& 0&2&3&2&2\\
       $\vert \vec \xi\vert^2 = 3$& 0&1&1&1&1\\
  \end{tabular}
  \caption[Multiplicities of $2\pi$ representations in the resting vector meson]{Multiplicities of $2\pi$ representations in the resting vector meson depending on the absolute momentum and taste. The multiplicities can be computed using the scalar product of characters of the representations of the vector states and the two-pion state, respectively.}
  \label{tab:multies}
\end{table}
\begin{table}
\centering
\begin{tabular}{|c||c|c|}
\hline
&$\vec \xi^2 = 0(3)$&$\vec \xi^2 = 1(2)$\\
\hline\hline
$\vec p^2 = 1$&$C^{1,\alpha}_{\lambda \vec e_i, \vec 0} = \lambda \frac{1}{\sqrt{2}}\delta^{\alpha i}$&$C^{1,\alpha}_{\lambda \vec e_i,\vec f_j} = \lambda \frac{1}{\sqrt{2}} \delta^{ij}\delta^{\alpha i}$\\
&&$C^{2,\alpha}_{\lambda \vec e_i,\vec f_j} = \frac{\lambda}{2}(1-\delta_{ij})\delta^{\alpha i}$\\
\hline
$\vec p^2 = 2$&$C^{1,\alpha}_{\lambda \vec e_i + \mu e_j,\vec 0} = \frac{\lambda}{2\sqrt{2}}\delta^{\alpha i}$&$C^{1,\alpha}_{\lambda \vec e_i + \mu \vec e_j,\vec f_k} = \frac{\lambda}{2\sqrt{2}}\vert\epsilon^{ijk}\vert\delta^{\alpha i}$ \\
&&$C^{2,\alpha}_{\lambda \vec e_i + \mu \vec e_j,\vec f_k} = \frac{\lambda}{2\sqrt{2}}\delta^{ki}\delta^{\alpha i}$ \\
&&$C^{3,\alpha}_{\lambda \vec e_i + \mu \vec e_j,\vec f_k} = \frac{\lambda}{2\sqrt{2}}\delta^{kj}\delta^{\alpha i}$ \\
\hline
$\vec p^2 = 3$&$C^{1,\alpha}_{\lambda \vec e_i + \mu \vec e_j + \nu \vec e_k,\vec 0} = \frac{\lambda}{2\sqrt{2}}\delta^{\alpha i}$&$C^{1,\alpha}_{\lambda \vec e_i + \mu \vec e_j + \nu \vec e_k,\vec f_l} = \frac{\lambda}{2\sqrt{2}}\delta^{il}\delta^{\alpha i}$ \\
&&$C^{2,\alpha}_{\lambda \vec e_i + \mu \vec e_j + \nu \vec e_k,\vec f_l} = \frac{\lambda}{4}\vert\epsilon^{jkl}\vert\delta^{\alpha i}$ \\
\hline
$\vec p^2 = 4$&$C^{1,\alpha}_{2\lambda \vec e_i, \vec 0} = \lambda \frac{1}{\sqrt{2}}\delta^{\alpha i}$&$C^{1,\alpha}_{2\lambda \vec e_i,\vec f_j} = \lambda \frac{1}{\sqrt{2}} \delta^{ij}\delta^{\alpha i}$\\
&&$C^{2,\alpha}_{2\lambda \vec e_i,\vec f_j} = \frac{\lambda}{2}(1-\delta_{ij})\delta^{\alpha i}$\\
\hline
\end{tabular}
\caption[Clebsch-Gordan coefficients for the vector states]{The Clebsch-Gordan coefficients needed for constructing the vector meson out of two-pion states.}
\label{tab:CG}
\end{table}

\section{Implementation}
\label{sec:implementation}
The coordinate-space correlation function of the $\pi^+$ operator is defined as
\begin{align}
  C_{\pi^+}(x,y) = [ u(x) O(x) \bar d(x) ][ d(y) O^{\dagger}(y) \bar u(y)],
\end{align}
whereas the $\pi^-$ correlation function is
\begin{align}
  C_{\pi^-}(x,y) = [ \bar u(x) O(x) d(x) ][ \bar d(y) O^{\dagger}(y) u(y)].
\end{align}
Here, we define $O(x) = \left(\gamma_5\otimes \Xi\right) \exp\left(-\mathrm{i}\vec p \cdot \vec x \right) $ in the spin-taste representation with a certain taste $\Xi$, the momentum is given by $\vec p = \frac{2\pi}{L}\vec n$ with $L$ being the spatial lattice extend and $\vec n \in \mathbb{N}^3$. The definition of the spin-taste operators for staggered fermions is given in Section~\ref{sec:stagg}. From now on I will drop the coordinate of the operator $O$ since it is given implicitly by the neighboring arguments.\\
To construct two-pion states we consider a combination of $\pi^+$ and $\pi^-$ with negative charge conjugation quantum number given by:
\begin{align}
  \frac{\pi^+\pi^- - \pi^-\pi^+}{2}
\end{align}
This is done in order to match the strong isospin and charge conjugation quantum number of the vector meson.\\
Therefore, the correlation function of the antisymmetrized $\pi^+\pi^-$-state has to be considered:
\begin{align*}
  C_{\pi^+\pi^-}(x,y,z,w) = [ u(x) O \bar d(x) ][ d(y) P \bar u(y) ][ d(z) O^{\dagger} \bar u(z)][ u(w) P^{\dagger} \bar d(w)]
\end{align*}
From this we get four different contractions:
\begin{align}
\label{eq:conn_contraction1}
 \langle  C_{\pi^+\pi^-}(x,y,z,w) \rangle =  &\underbrace{\textrm{Tr}\left[M^{-1}_u(y,x)OM^{-1}_d(x,y)P\right] \textrm{Tr} \left[ M^{-1}_d(w,z)O^{\dagger}M^{-1}_u(z,w)P^{\dagger}\right]}_{\textrm{disconnected, ``t-channel''}}\\
  \label{eq:conn_contraction2}
  -& \underbrace{\textrm{Tr}\left[M^{-1}_u(z,x)OM^{-1}_d(x,y)PM^{-1}_d(w,z)O^{\dagger}M^{-1}_u(y,w)P^{\dagger}\right]}_{\textrm{connected}}\\
  \label{eq:conn_contraction3}
  -& \underbrace{\textrm{Tr}\left[M^{-1}_u(y,x)OM^{-1}_d(x,z)PM^{-1}_d(w,y)O^{\dagger}M^{-1}_u(z,w)P^{\dagger}\right]}_{\textrm{connected}}\\
  \label{eq:conn_contraction4}
   +& \underbrace{\textrm{Tr}\left[M^{-1}_u(z,x)OM^{-1}_d(x,z)P\right] \textrm{Tr} \left[M^{-1}_d(w,y)O^{\dagger}M^{-1}_u(y,w)P^{\dagger}\right]}_{\textrm{disconnected, ``s-channel''}}
\end{align}
Here, $\langle ...\rangle$ denotes the expectation value with respect to the fermion path-integral at a fixed gauge configuration. The trace is taken with respect to the color components. As mentioned before, it is necessary to antisymmetrize the ingoing as well as the outgoing pions. So we are left with 
\begin{align}
 \frac{1}{4}\left( C_{\pi^+\pi^-}(x,y,z,w) - C_{\pi^+\pi^-}(y,x,z,w) - C_{\pi^+\pi^-}(x,y,w,z) - C_{\pi^+\pi^-}(y,x,w,z)\right).
\end{align}
The first disconnected term (Equation~\eqref{eq:conn_contraction1}) will be zero due to its symmetric structure in combination with the antisymmetrization of the two-pion correlator. Therefore, no implementation for this term is needed and the term in Equation~\eqref{eq:conn_contraction4} will be called \textit{disconnected} from now on. Additionally, one will see that Equation~\eqref{eq:conn_contraction2} and Equation~\eqref{eq:conn_contraction3} are indistinguishable after the antisymmetrization. Therefore, we will only consider one of those and call it the \textit{connected} diagram. One has to note that this connected diagram then gets a prefactor of $-2$. Within the analysis one should ensure that a proper antisymmetrization with respect to the ingoing and outgoing pions is done.\\
Since the pions need to have back-to-back momentum one has to add certain phase factors before summing up the spatial components. We consider the center-of-mass frame and therefore, the ingoing and the outgoing momenta have to sum up to zero. The momentum of each ingoing pion is $\vec p$ (resp. $-\vec p$) and the momentum for the outgoing pions is $\vec q$ (resp. $-\vec q$).\\
Of course, the two-pion states are not the only states needed in the GEVP for the reconstructions of the vector meson. One also has to implement the vector meson ($\rho$-meson) correlation function it self to resolve all energy eigenstates. This correlation function is given by
\begin{align}
  C_{\rho}(x,z) &= \frac{1}{2}\left([ u(x) O \bar u(x) ] - [ d(x) O^{\dagger} \bar d(x)]\right)\left([ \bar u(z) O u(z) ] - [  \bar d(z) O^{\dagger} d(z)]\right),
\end{align}
where $O = \gamma_i \otimes \mathbbm{1}$ in the spin-taste representation, which corresponds to a taste-singlet vector meson. Finally, one has only two contributing contractions which corresponds to one contributing contraction in the isospin-symmetric limit, which is used for the implementation:
\begin{align}
\Tr\left[ M^{-1}_{u/d}(x,z)OM^{-1}_{u/d}(z,x)O^{\dagger}\right] 
\end{align}
The GEVP also requires terms that mix the vector and the two-pion state of which the correlation functions are given by:
\begin{align}
   C_{\rho\pi\pi}(x,z,w) &= \frac{1}{2}\left([ u(x) O \bar u(x) ] - [ d(x) O^{\dagger} \bar d(x)]\right)[ d(z) P^{\dagger} \bar u(z)][ u(w) P^{\dagger} \bar d(w)].
\end{align}
In comparison to the pure pion or vector states, this three-point correlation functions are not real but complex valued. In addition one can prove the identity
\begin{align}
   C_{\rho\pi\pi}(x,w,z) = C_{\rho\pi\pi}(x,z,w)^*,
\end{align}
which, combined with the antisymmetrization in the pions, leads to the following expression that has to be implemented:
\begin{align}
  \Im \Tr \left[ M^{-1}_{u/d}(w,z)OM^{-1}_{u/d}(x,w)P^{\dagger}M^{-1}_{d/u}(z,w)P\right].
\end{align}
For the implementation, we use certain noise sources $\xi_{e/o,t_0}$ that fulfill the following properties:
\begin{align}
  & \xi_{e,t_0}^T\xi_{e,t_1} = \frac{N_c L^3}{2}\delta_{t_0,t_1},\\ 
  & \xi_{e,t_0}^T\xi_{e,t_1} = 0,\\ 
  & \xi_{e,t_0}^T\xi_{e,t_1} = 0,\\ 
  & \xi_{e,t_0}^T\xi_{e,t_1} = \frac{N_c L^3}{2}\delta_{t_0,t_1}.
\end{align}
In this case, they are represented by wall sources in the following way:
\begin{align}
\label{eq:wallsource}
\begin{aligned}
    \xi_{e,t_0}(\vec x,t,c) = \begin{cases}
    0 \quad &\textrm{if $\vert x \vert_1 \%2 = 1$}\\
    (-1)^{r(\vec x,c)} \delta_{t_0,t}\quad &\textrm{else}\\
  \end{cases},\\
    \xi_{o,t_0}(\vec x,t,c) = \begin{cases}
    0 \quad &\textrm{if $\vert x \vert_1 \%2 = 0$}\\
    (-1)^{r(\vec x,c)}\delta_{t_0,t} \quad &\textrm{else}\\
  \end{cases}. 
\end{aligned}
\end{align}
Here, $r(\vec x,c)$ is a randomly chosen integer for each spatial lattice site and each color index. \\
For certain cases (vector mesons) we need point sources instead of the aforementioned wall sources. Those can be written as
\begin{align}
\label{eq:pointsource}
  \xi_{e/o,t_0}(\vec x,t,c) = \sqrt{\frac{N_cL^3}{2}}\delta(\vec y - \vec x)\delta_{t_0,t}\delta_{d,c}
\end{align}
Here $\vec y$ is a randomly chosen spatial lattice site (that is even ($e$) or odd($o$)) and $d$ is a random color index.
\begin{figure}[h]
\begin{center}
\begin{tikzpicture}[very thick,q0/.style={->,DarkBlue,semithick,yshift=5pt,shorten >=5pt,shorten <=5pt}]

 % Loop
 \def\radius{0.9}
 
 \draw[] (0,-\radius) -- (2*\radius,-\radius) ;
 \draw[] (0,-\radius) -- (0,\radius) ;
 \draw[] (2*\radius,-\radius) -- (2*\radius,\radius) ;
 \draw[] (0,\radius) -- (2*\radius,\radius) ;
  \node[left] (1) at (0,0) {$u$};
    \node[right] (1) at (2*\radius,0) {$u$};
  \node[above] (1) at (\radius,\radius) {$\bar d$};
  \node[below] (1) at (\radius,-\radius) {$\bar d$};
  \node[above] (1) at (-0.2*\radius,\radius) {$\pi^{\dagger}(\vec x,0)$};
  \node[below] (1) at (-0.2*\radius,-\radius) {$\pi(\vec y,0)$};
  \node[above] (1) at (2.2*\radius,\radius) {$\pi(\vec z,t)$};
  \node[below] (1) at (2.2*\radius,-\radius) {$\pi^{\dagger}(\vec w,t)$};
  \node[below] (1) at (\radius,-1.8*\radius) {$a)$};
\def \offset{3.5}
\draw[] (\offset,-\radius) -- (\offset+2*\radius,-\radius) ;
 \draw[] (\offset,-0.9*\radius) -- (\offset+2*\radius,-0.9*\radius) ;
 \draw[] (\offset,\radius) -- (\offset+2*\radius,\radius) ;
 \draw[] (\offset,0.9*\radius) -- (\offset+2*\radius,0.9*\radius) ;
  \node[below] (1) at (\offset+\radius,\radius) {$u$};
  \node[above] (1) at (\offset+\radius,\radius) {$\bar d$};
  \node[below] (1) at (\offset+\radius,-\radius) {$u$};
  \node[above] (1) at (\offset+\radius,-\radius) {$\bar d$};
  \node[above] (1) at (\offset-0.2*\radius,\radius) {$\pi^{\dagger}(\vec x,0)$};
  \node[below] (1) at (\offset-0.2*\radius,-\radius) {$\pi(\vec y,0)$};
  \node[above] (1) at (\offset+2.2*\radius,\radius) {$\pi(\vec z,t)$};
  \node[below] (1) at (\offset+2.2*\radius,-\radius) {$\pi^{\dagger}(\vec w,t)$};
    \node[below] (1) at (\offset + \radius,-1.8*\radius) {$b)$};

\begin{feynman}
  \vertex (a) at (\offset+0.6*\radius,\radius);
  \vertex (b) at (\offset+1.7*\radius,-\radius);
  \diagram*{
   (a) -- [gluon] (b);
  };
\end{feynman}



  \draw[] (2*\offset,-\radius) -- (2*\offset,+\radius) ;
 \draw[] (2*\offset + 0.1*\radius,-\radius) -- (2*\offset + 0.1*\radius, +\radius) ;
 \draw[] (2*\offset+2*\radius,-\radius) -- (2*\offset+2*\radius,+\radius) ;
 \draw[] (2*\offset + 1.9*\radius,-\radius) -- (2*\offset + 1.9*\radius, +\radius) ;
  \node[left] (1) at (2*\offset,0) {$u$};
  \node[right] (1) at (2*\offset,0) {$\bar d$};
  \node[left] (1) at (2*\offset+2*\radius,0) {$u$};
  \node[right] (1) at (2*\offset+2*\radius,-0) {$\bar d$};
 \begin{feynman}
  \vertex (a) at (2*\offset,0.6*\radius);
  \vertex (b) at (2*\offset+2*\radius,-0.6*\radius);
  \diagram*{
   (a) -- [gluon] (b);
  };
  \end{feynman}
  
  \node[above] (1) at (2*\offset-0.2*\radius,\radius) {$\pi^{\dagger}(\vec x,0)$};
  \node[below] (1) at (2*\offset-0.2*\radius,-\radius) {$\pi(\vec y,0)$};
  \node[above] (1) at (2*\offset+2.2*\radius,\radius) {$\pi(\vec z,t)$};
  \node[below] (1) at (2*\offset+2.2*\radius,-\radius) {$\pi^{\dagger}(\vec w,t)$};
      \node[below] (1) at (2*\offset + \radius,-1.8*\radius) {$c)$};

\end{tikzpicture}

\begin{tikzpicture}[very thick,q0/.style={->,DarkBlue,semithick,yshift=5pt,shorten >=5pt,shorten <=5pt}]

 % Loop
 \def\radius{0.9}
 
 \draw[] (0,0) -- (2*\radius,-\radius) ;
 \draw[] (2*\radius,-\radius) -- (2*\radius,\radius) ;
 \draw[] (0,0) -- (2*\radius,\radius) ;
  \node[right] (1) at (2*\radius,0) {$\bar d$};
  \node[above] (1) at (\radius,0.5*\radius) {$u$};
  \node[below] (1) at (\radius,-0.5*\radius) {$u$};
  \node[left] (1) at (0,0) {$\rho^{\dagger}(\vec x,0)$};
  \node[above] (1) at (2.2*\radius,\radius) {$\pi(\vec z,t)$};
  \node[below] (1) at (2.2*\radius,-\radius) {$\pi^{\dagger}(\vec w,t)$};
      \node[below] (1) at (\radius,-1.8*\radius) {$d)$};
  \def\offset{5*\radius}
 \draw[] (\offset,0) -- (\offset + 2*\radius,0) ;
 \draw[] (\offset,0.1) -- (\offset + 2*\radius,0.1) ;
 \draw[] (0,0) -- (2*\radius,\radius) ;
  \node[above] (1) at (\offset +\radius,0.2*\radius) {$u$};
  \node[below] (1) at (\offset +\radius,-0.2*\radius) {$\bar d$};
  \node[left] (1) at (\offset +0,0) {$\rho^{\dagger}(\vec x,0)$};
  \node[left] (1) at (\offset +3.5*\radius,0) {$\rho(\vec z,0)$};
    \node[below] (1) at (\offset + \radius,-1.8*\radius) {$e)$};

\end{tikzpicture}
\end{center}
\caption[Diagrams entering the vector-vector GEVP]{The diagrams corresponding to the GEVP including vector-meson-type states. The first one (conn., $a)$) includes interchange of quarks in the two pion states, the second one (free, $b)$) describes two pions propagating independently and the third one (disc., $c)$) is identical to zero due to the isospin combination of the pions. The decay of the vector mesons into two pions is described by the triangle ($d)$) and the remaining diagram ($e)$) is the vector meson correlation function~\cite{Frech:2023bzy}}
\label{fig:diagrams}
\end{figure}
\subsection{Vector meson correlator}
\label{subse:vectormesoncorrelator}
% The vector meson correlator has to be included with unsmeared as well as smeared operators. Therefore, it is sufficient to use a point source instead of wall source. For the quark smearing we use Wuppertal smearing as it is done in~\cite{Borsanyi_2021}. The action Wuppertal smearing operator on a lattice vector $v$ is given by
% \begin{align}
%   Wv(x) = (1-\sigma)v(x) +\\
%   \frac{\sigma}{6}\sum_{\mu = 1}^3\left (U_{\mu}(x)U_{\mu}(x + \hat \mu)v(x + 2\hat \mu)+ U^{\dagger}_{\mu}(x-\hat\mu)U^{\dagger}_{\mu}(x -2 \hat \mu)v(x -
%    2\hat \mu)\right)
% \end{align}
% where $\sigma$ is the smearing coefficients. The gauge links included in this operator are also smeared before using a three dimensional stout smearing procedure~\cite{Morningstar_2004}~\cite{Borsanyi_2021}:
% \begin{align}
%   \tilde U_i(x) = \exp\left(iQ_{i}(x) \right)U_i(x)\\
%   Q_i(x) = \frac{i}{2}\left(\Omega_i^{\dagger}(x) - \Omega_i(x) \right) - \frac{i}{6}\Tr\left(\Omega_i^{\dagger}(x) - \Omega_i(x) \right)\\
%   \Omega_i(x) = C_i(x)U_i^{\dagger}(x)\\
%   C_i(x) = \rho \sum_{i \neq j}\left(U_j(x)U_i(x + \hat j)U_j^{\dagger}(x+\hat i) + U_j^{\dagger}(x-\hat j)U_i(x-\hat j)U_j(x-\hat j+\hat i) \right)
% \end{align}
% The smearing coefficients and the number of application of the smearing operator in this study are given in table \ref{tab:smearing_coeff}.
% \begin{table}[]
%   \centering
%   \begin{tabular}{c|c|c|c|c|c|c}
%     $\beta$ & $L/a$ & $T/a$ & $\sigma$ & $N_{\textrm{Wuppertal}}$&$ \rho$&$N_{\textrm{stout}}$ \\
%     \hline
%      3.7000&64&32& 0.125&12&0.125&32\\
%   \end{tabular}
%   \caption{Smearing coefficients for the Wuppertal smearing on the vector meson}
%   \label{tab:smearing_coeff}
% \end{table}

The computation of the vector meson correlation function is done by implementing the scalar product
\begin{align}
\label{eq:rho_corr}
  C_{\rho}(t_0,t) = \frac{1}{N_{\xi}}\sum_{\xi}\sum_{t_1,t_2}\braket{\xi_{t_0} M^{-1}_{t_0,t_1} B \left(P_{t}\right)_{t_1,t_2}}{ M^{-1}_{t_2,t_0} A\xi_{t_0}} ,
\end{align}
where $A/B$ denote components of the taste-singlet vector meson operator (in this case we have $A = B$) and $P_t$ is a projection on a fixed time slice $t$ that is given by
\begin{align}
  \left(P(t)\right)_{t_1,t_2} = \delta_{t,t_1}\delta_{t,t_2}.
\end{align}
From now on, this projection is included implicitly. The expectation value is computed by the mean values with different random sources $\xi$ that are only non-zero at a randomly chosen time slice $t_0$ according to Equation \eqref{eq:wallsource} and \eqref{eq:pointsource}. The number of inversion, which is the most costly factor of such a computation for the vector meson correlator is given by
\begin{align}
  N_\xi(1+N_O) = N_\xi(1+3)= 4N_\xi
\end{align}
where $N_\xi$ gives the number of random sources $\xi$ and $N_O$ is the number of vector meson operators (smeared and unsmeared, three components). Unfortunately, this expression converges only very slowly towards the exact result when increasing $N_\xi$. Therefore, for our usual values of $N_\xi$ the vector meson correlator is dominated by noise. In order to get rid of the problem, we include \textit{low-mode averaging} (LMA). The idea behind LMA is that we expect the correlator mainly to depend on the eigenvalues of $M$ with small absolute value. From the preconditioning, we already know the lowest eigenvalues $\lambda_i$ and the corresponding eigenvectors $v_i$ so we can compute the most contributing part exactly:
\begin{align}
 C^{LMA}_{\rho}(t_0,t) = \sum_{ij}^{N_\text{eig}}\frac{1}{\lambda_i\lambda_j} \bra{v_i O_A} P_{t_0}\ket{ v_j}\bra{v_j O_B }P_t\ket{ v_i}.
\end{align}
The eigen-part of the inversion is defined as
\begin{align}
  M_e^{-1} = \sum_i^{N_\text{eig}} \frac{1}{\lambda_i}\ketbra{v_i}{v_i}
\end{align}
and is of course, not the true inversion of the staggered matrix. There is a rest-part we define as
\begin{align*}
  M_r^{-1} = M^{-1} - M_e^{-1}.
\end{align*}
Due to the fact that the correlator in Equation \eqref{eq:rho_corr} includes two inversions of the staggered matrix, we have eigen-eigen (ee), rest-rest (rr) and eigen-rest and rest-eigen (er) contributions in the vector meson correlator. If the number of pre-calculated eigenvectors gets to a suitable extend, the eigen-eigen part is the dominant one. This part then will be replaced by the low-mode averaging expression to receive a strong improvement of the signal-to-noise ratio especially for larger Euclidean time distances. The eigen-eigen part of Equation \eqref{eq:rho_corr} can be used to cross-check the low-mode-averaging. 
\subsection{Free pion diagram}
\label{sec:free_pion}
The free pion diagram has to advantage that it can be computed by multiplying two one-pion correlation functions. The computation looks similar as the vector meson in the previous subsection, but the operator has to be replaced by the corresponding pion-operator with proper taste quantum numbers. Due to the fact that the pion signal is very clean (since we deal with $\pi^+/\pi^-$) low-mode averaging is not necessary. In addition we do not need to smear the operators and, therefore, we can use a wall-source instead of a point source. The main problem in the implementation of this diagram is that we might have time-non-local operators. In the lattice simulation, an operator includes shift by $(s_{\mu} + t_{\mu})\mod 2$ in the $\mu$ direction. Whenever $\mu$ is spatial this is no problem, because the shift in symmetrized in forward and backward direction. But since we measure correlation functions at a given time slice, we cannot do this in time direction. The application of a time-non-local operator would shift the time slice at which the operator is measured from $t$ to $t-1$. Therefore, we also include sources at $t_0 + 1$ and average over different correlation functions in the following way, where $O^l$ denotes a local and $O^n$ a time-non-local operator. We start explaining the case in which $A$ is a local operator and $B$ can be either local or non-local. The first case is well-known:
\begin{align}
   C_{B^l A^l}(t_0,t) \coloneqq \langle \xi_{t_0} M^{-1}_{t_0,t} B^l \vert M^{-1}_{t,t_0} A^l\xi_{t_0} \rangle,\\
   C_{B^n A^l}(t_0,t) \coloneqq \frac{1}{2}\left( \langle \xi_{t_0} M^{-1}_{t_0,t} B^n \vert M^{-1}_{t,t_0} A^l \xi_{t_0}\rangle + \langle \xi_{t_0} M^{-1}_{t_0,t} \vert \left(B^n\right)^{\dagger} M^{-1}_{t,t0} A^l\xi_{t_0} \rangle\right).
\end{align}
At the vertical line we compute the scalar-product which includes a projection on time slice $t$, but the difference of $t$ and $t_0$ differs by one between the two terms in $  C_{B^n A^l} $. To include non-local operators for $A$, we compute
\begin{align*}
   C_{B^l A^n}(t_0,t) \coloneqq \frac{1}{2}\left( \langle \xi_{t_0} M^{-1}_{t_0,t} B^l \vert M^{-1}_{t,t_0} A^l\xi_{t_0} \rangle + \langle \xi_{t_0+1}\left( A^n\right)^{\dagger} M^{-1}_{t_0,t} B^l \vert M^{-1}_{t,t_0+1} \xi_{t_0+1} \rangle\right)\\
\end{align*}
Combining the two methods for computing non-local $A$ resp. $B$ one yields the expression for both to be non-local:
\begin{align}
\begin{aligned}
    C_{B^n A^n}(t_0,t) \coloneqq \frac{1}{4}\bigg( \langle \xi_{t_0} M^{-1}_{t_0,t} B^l \vert M^{-1}_{t,t_0} A^l\xi_{t_0} \rangle + \langle \xi_{t_0+1}\left( A^n\right)^{\dagger} M^{-1}_{t_0+1,t} B^l \vert M^{-1}_{t,t_0+1} \xi_{t_0+1} \rangle &\\
  + \langle \xi_{t_0} M^{-1}_{t_0,t} B^n \vert M^{-1}_{t,t_0} A^l \xi_{t_0}\rangle + \langle \xi_{t_0} M^{-1}_{t_0,t} \vert \left(B^n\right)^{\dagger} M^{-1}_{t,t_0} A^l\xi_{t_0} \rangle + \langle \xi_{t_0 + 1} M^{-1}_{t_0 + 1,t} B^l \vert M^{-1}_{t,t_0+1} A^l\xi_{t_0+1} \rangle &\bigg).
\end{aligned}
\end{align}
The free part will be the dominating contribution to the two-pion interaction especially for small time separations. Nevertheless for larger separations, the connected part may play a more significant role.
\subsection{Interaction diagram}
For the implementation of the interaction diagram between the $\rho$-meson and the two-pion state, there are two possible ways of implementation. Let $A$ be the operator describing one of the corresponding pions. It has to be mentioned that both $A$ operators will have opposite momentum, so one has to take care of this issue when implementing the diagrams. The $\rho$ operator will be denoted by $B$. The first method is the implementation of the decay:
\begin{align}
  C_{\rho A A}(t_0,t)  = \langle\xi_{t_0}M^{-1}_{t_0,t}A \textcolor{red}{M^{-1}_{t,t}}A\vert M^{-1}_{t_0,t_0}\rho\xi_{t_0}\rangle
\end{align}
Here, a problem appears that also arises in the next section that discusses the connected pion diagram. The red inversion has to be perform on each time-slice which will increases the computational time drastically. Luckily, one can compute the diagram in the opposite direction:
\begin{align}
    C_{AA\rho} = \langle\xi_{t_0}M^{-1}_{t_0,t}B\vert M^{-1}_{t,t_0}AM^{-1}_{t_0,t_0}A\xi_{t_0}\rangle.
\end{align}
This way, the dominant part of computational costs comes from the right-hand-side of the scalar product where it is proportional to the number of operators included but not to the number of time-slices.
\subsection{Connected pion diagram}
The connected part of the diagram is by far the most expensive part of the computation and requires the most elaborated strategy for a efficient computation. The main issue for the high costs can be seen from a naive implementation for local operators:
\begin{align}
  C^\text{conn.}_{B B A A}(t_0,t)  = \langle\xi_{t_0}M^{-1}_{t_0,t}B \textcolor{red}{M^{-1}_{t,t}}B\vert M^{-1}_{t,t_0}AM^{-1}_{t_0,t_0}A\xi_{t_0}\rangle .
\end{align}
The red inversion has to be performed on every timeslice, so that there occurs an overall number of inversions proportional to~\cite{Frech:2025fjp}
\begin{align}
  T/a \times \#O \times \# N_\xi
\end{align}
. This is no problem on small volumes but since the cost of a single inversion scales roughly with $V = L^3 \times T$, the additional time scaling becomes a problem at larger volumes. An approach to save computational time for the connected diagrams was presented by me at the LATTICE2024 conference. The key ingredient is the low-mode averaging that was also used in Section \ref{subse:vectormesoncorrelator}. Therefore, the red inversion in replaced by the sum of the eigen-part $M_e^{-1} = \sum_i^{N_\text{eig}} \frac{1}{\lambda_i}\ket{v_i}\bra{v_i}$ and $M_r^{-1} = M^{-1} - M_e^{-1}$:
\begin{align}
  C^\text{conn.}_{B B A A} &= \langle\xi_{t_0}M^{-1}_{t_0,t}B \sum_i^{N_\text{eig}} \frac{1}{\lambda_i}\ket{v_i}\bra{v_i}B\vert M^{-1}_{t,t_0}AM^{-1}_{t_0,t_0}A\xi_{t_0}\rangle\\
  &+ \langle\xi_{t_0}M^{-1}_{t_0,t}B \left(M^{-1}_r\right)_{t,t}B\vert M^{-1}_{t,t_0}AM^{-1}_{t_0,t_0}A\xi_{t_0}\rangle.
\end{align}
Of course the rest part also requires inversions at each time slice. Therefore, an approximation of the identity $\mathbbm{1} \approx \sum_\sigma \ket{\sigma} \bra{\sigma}$ is included, where the vectors $\ket{\sigma}$ are normalized stochastic sources.
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/compareplot.pdf}
    \caption[Comparison of the quality, exact inversion vs. low-mode-averaging]{Comparison of exact and estimated two-pion propagator on the ensemble with $\beta = 3.7000$, $L/a = 48$ and $T/a = 64$. Left panel: Correlation function of the two-pion propagator with both pions having Goldstone and two units of momentum (in opposite directions). The blue data are measured exactly using the full inversion whereas the orange are exact on the lowest 1500 eigenvalues with the rest part being estimated using 24 random sources. Right panel: Ratio of the propagators from the two different methods. For a better visibility the optimal ratio is marked by a blue line. Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:compare_quality}
\end{figure}
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/scalingplot.pdf}
    \caption[Comparison of the computational time, exact computation vs low-mode-averaging]{Comparison of the computational time for measuring all observables on a single configuration using the full inversion and the low-mode-averaging, respectively. Left panel: Computational time depending on the volume. Right panel: Computational time depending on the number of operators. Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:compare_time}
\end{figure}
The number of stochastic sources $N_\sigma$ has to be fine-tuned, not least because we are interested in saving computational resources. Of course, we also need a sufficient large amount of random sources such that we can reach the desired accuracy. It could be shown that choosing $N_\sigma \sim L/(2a)$ is an appropriate first choice.\\
The quality of the result and the savings of computational time are the two things that need further investigations. In order to do this we employ representative configurations out of the repository for the $g-2$ project, a list of the ensembles can be found in Table~\ref{tab:configs} and Table~\ref{tab:configs_fs}. Finally, it was made the decision to explore the coarsest ensembles available ($\beta = 3.7000$). This has the advantage that we can save computational resources even if we do simulations at our usual volume of $(6^3\times 9)\,\mathrm{fm}^4$. Additionally, we have access to different volumes at this $\beta$ and therefore are able to investigate the volume scaling of the computing time.\\
When comparing the quality of the exact method and the approximation with low-mode-averaging and random sources it could be seen, that the quality of the second method strongly depends on the operator. For the  Goldstone taste in combination with zero-momentum it is fairly easy to get a good quality even with a really small choice for $N_\sigma$. This is related to the fact that in this case the two-pion propagator is dominated by the low-modes. However, $p = 0$ does not have an overlap to the vector states because it has no particular direction. Thus, one cannot use this method in the way it works best.\\ 
Unfortunately, the quality gets worse and worse the higher the momentum and the more non-local the taste. This is due to the fact that in these cases the contribution of higher modes starts to dominate the overall result. Increasing the number of eigenvalues and the number of random sources are possible solutions, but of course they come along with more computational effort. As a representative example I show a comparison of the quality of the Goldstone taste with the momentum being $\frac{4\pi}{L}$, i.e., two times the smallest possible momentum. It is shown in Figure~\ref{fig:compare_quality}. The propagators in the right panel look very similar apart from the noisy region around $T/2$. The comparison in the right panel is much more useful in this case. Here, the ratio of exact and low-mode-averaged is given. It is evident that the deviation starts to be significant around $18 a \approx 2.3\,\mathrm{fm}$. However, this is not a big deal since we want to use early windows to extract the mass plateaus and matrix elements. Thus larger noise at later windows does not provide any problems.\\
Apart from the fact that the new method seems to be sufficiently precise, the question remains whether we can save computational time. In order to clarify this, two tests were performed whose results are given in Figure~\ref{fig:compare_time}. The first measures the simulation time as a function of the volume. In both cases, a behavior of polynomial order can be observed with the low-mode-averaging method scaling slightly better. The other test is performed with respect to the number of operators included. Here, it can be seen that both methods also scale similarly, with the low-mode-averaged method being superior by a global factor of three. All these tests are performed with the number of eigenvalues and $N_\sigma$ being tuned so that both methods give a similar precision on the result.


\section{Decomposition of the vector-like states and reconstruction of the vector propagator}
\label{sec:tail_recon}
In this section, a few results on the tail reconstruction of the vector meson propagator are discussed. We use 48 configurations of the $\beta = 3.7000$ ensemble with a physical extension of $(6^3\times 9)\,\mathrm{fm}^4$. The lattice spacing is (depending on the scheme) slightly larger than $0.13\,\mathrm{fm}$ and we use $2+1+1$ quarks that are close to the (isospin-symmetric) physical point. The configurations are generated with the \texttt{4stout} action~\cite{Borsanyi_2021,Boccaletti:2024guq,FKAPI,Morningstar_2004,Luscher:1984xn}.\\
On this configuration, the two-pion propagators (connected and disconnected), the vector correlator and the interaction propagator are measured. Of course, not all possible two-pion states are included in this study, but only those with an energy smaller than the mass of the vector meson $m_\rho$. This is the case for Goldstone pions with up to three units of momentum, for axial-taste pions with up to two units of momentum, and for tensor-taste pions with one unit of momentum. For nearly all of the relevant combinations, 120 random sources are used. However, since the connected diagram is the most expensive and, additionally, very small compared to the disconnected one, only 30 random sources are used.\\


To resolve the energy states that contribute to the vector-meson propagator, a generalized eigenvalue problem is constructed. To deal with the noise related to higher excitation of the vector meson ($E > m_\rho$), we use the pencil-of-functions method~\cite{DeTar_2015,370583}. This means that time-shifted $\rho$-propagators are also included in the GEVP. This method is also applied in Chapter~\ref{cha:LatticeTechniques}. As it can be seen in Table~\ref{tab:CG} there are now ten different two-pion combinations with an energy lower than that of the $\rho$-meson. Together with tree different time-shifts of the vector propagator from the pencil-of-functions method the overall correlation matrix is given by
\begin{align}
    \mathbf{C}(t) = \begin{pmatrix}
        C_{\rho\rho}(t-4)&C_{\rho\rho}(t-3)&C_{\rho\rho}(t-2)&C_{\rho(\pi\pi)^1}(t-2)&C_{\rho(\pi\pi)^2}(t-2)&...\\
        C_{\rho\rho}(t-3)&C_{\rho\rho}(t-2)&C_{\rho\rho}(t-1)&C_{\rho(\pi\pi)^1}(t-1)&C_{\rho(\pi\pi)^2}(t-1)&...\\
        C_{\rho\rho}(t-2)&C_{\rho\rho}(t-1)&C_{\rho\rho}(t)&C_{\rho(\pi\pi)^1}(t)&C_{\rho(\pi\pi)^2}(t)&...\\
        C_{(\pi\pi)^1\rho}(t-2)&C_{(\pi\pi)^1\rho}(t-1)&C_{(\pi\pi)^1\rho}(t)&C_{(\pi\pi)^1(\pi\pi)^1}(t)&C_{(\pi\pi)^1(\pi\pi)^2}(t)&...\\
         C_{(\pi\pi)^2\rho}(t-2)&C_{(\pi\pi)^2\rho}(t-1)&C_{(\pi\pi)^2\rho}(t)&C_{(\pi\pi)^2(\pi\pi)^1}(t)&C_{(\pi\pi)^2(\pi\pi)^2}(t)&...\\
         ...&...&...&...&...&...&
    \end{pmatrix}.
\end{align}
Here, we denote the two-pion states that have the same quantum numbers as the $\rho$-meson by $(\pi\pi)^i$, where we used $1 \leq i \leq 10$ to label the states whose total energy is bound by $770\,\mathrm{MeV}$.\\ 
We solve the GEVP
\begin{align}
     \mathbf{C}(t)\cdot V(t,t_0) =  \lambda(t,t_0)\mathbf{C}(t)\cdot V(t,t_0),
\end{align}
where we choose $t_0 = 5$ and $t = 6$. The correlators are reconstructed by \begin{align}
    \tilde \lambda_i(t) = V_i^\dagger\cdot \mathbf{C}(t) \cdot V_i.
\end{align}
They are shown in Figure~\ref{fig:vec_corr}. It is evident that there are two states that are very noisy, namely the green and the orange one. They are related to the higher excitations that appear in the $\rho$-propagator and are filtered out using the pencil of functions. Those two states are not considered in the reconstruction. However, this is not necessary because higher-energy states dominate the short-distance contributions and are suppressed exponentially in the long range. To get the energy states that contribute to the long distance contributions, we compute the corresponding effective masses using the local mass extraction formula
\begin{align}
\label{eq:effmasslocal2}
 M^{\textrm{eff, loc}}(t)= \frac{1}{\Delta}\cosh^{-1}\frac{G(t + \Delta) + G(t-\Delta)}{2G(t)},
 \end{align}
 where $\Delta$ is chosen to be $2a$.
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/corr.pdf}
    \caption[Eigenvalues of the vector GEVP]{Eigenvalues of the vector GEVP. Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:vec_corr}
\end{figure}
The corresponding effective masses are shown in Figure~\ref{fig:vec_meff}. Here, also the effective masses from the higher states (orange and green) are shown and it can be seen that they do not form a proper plateau at this level of precision. However, it is very likely that these states indeed correspond to higher energies, since most of the orange and green points are above the threshold of $am_\rho\approx 0.5$. 
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/E_eff.pdf}
    \caption[Effective masses of the vector GEVP]{Effective masses of the vector GEVP. Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:vec_meff}
\end{figure}
To reconstruct the matrix elements that are needed to reconstruct the vector meson correlation function, we use the expression
\begin{align}
    a_i^\text{eff.} = \frac{(V_i \cdot\mathbf{C}(t) \cdot e_\rho)^2}{\tilde \lambda_i(t)}\times \exp (m^\text{eff}(t) \times t),
\end{align}
where $e_\rho$ is the unit vector that picks out the unshifted $\rho$-components of the of the correlation matrix $\mathbf{C}(t)$.
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/states.pdf}
    \caption[Matrix elements of the GEVP]{Matrix elements of the GEVP.  Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:vec_ai}
\end{figure}
The corresponding effective matrix elements can be seen in Figure~\ref{fig:vec_ai}. Again, the green and the orange do not build a plateau; however, they are not considered in the reconstruction. Furthermore, apart from the blue state, all matrix elements have a very large relative error. Additionally, the blue state dominates the overall result by far. This suggests that there appeared some problems in the generalized eigenvalue problem. However, we want to use the data to reconstruct the first moment $\sum_t\frac{t^4}{24}\rho(t)$. The integrand of the first moment can be seen in Figure~\ref{fig:vec_rec}, where the reconstruction is compared to the result of the direct simulation of the vector propagator. The two data sets are not compatible in the short-distance window, however, this is due to the fact that we did not include the (very noisy) higher excitations. Nevertheless, in the longer distances, the two methods are compatible with the reconstruction being much more precise. 
\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/6/reconstruction.pdf}
    \caption[Reconstruction of the integrand of the first moment]{Reconstruction of the integrand of the first moment. Previously published in~\cite{Frech:2025fjp}.}
    \label{fig:vec_rec}
\end{figure}
\section{Summary and comment on the further progress}
\label{sec:conclusion}
In this work, a mathematical framework to to construct the two-pion states with the same quantum numbers as the taste-singlet vector meson at rest was derived. In principle, this framework can also be used for other tastes and momenta of the vector meson. The diagrams corresponding to this project (vector-meson, two-pion (connected and disconnected) and interaction) were implemented in the software framework \textit{dynqcd}. For the connected two-pion contributions, a method for more efficient computations was designed. Some tests were performed to compare this method to the direct computation of the connected diagram and it was found out that the new method is superior. Finally, all this methods were combined to reconstruct the vector meson on our coarsest ensemble with a size of $(6^3\times 9)\,\mathrm{fm}^4$, a lattice spacing of $a \approx 0.134\,\textrm{fm}$ and physical quark masses.\\
A precise computation of the tail of the vector-meson propagator is a very challenging task. Especially if staggered fermions are used, a lot of problems can occur. Staggered fermions have, in contrast to e.g. Wilson fermions, a taste structure. Thus, a large variety of states is present at finite lattice spacing, which have to be included in the generalized eigenvalue problem. For example, in this study ten two-pion states are used in the GEVP, without the taste structure it would have been three. Apart from the large number of contributing states, the taste structure provides more problems: The non-local tastes lead to an increment of the noise of the corresponding correlators, which makes it harder to resolve the energy states properly. As discussed in Section~\ref{sec:implementation}, the large number of operators require additional methods like low-mode averaging. It was also shown that this method suffer from non-locality of the tastes.\\
In addition to the large number of problems with the taste structure, solving the GEVP does not seem to work optimally in this case. The overlap of the $\rho$ and the two-pion state is likely too small. This issue could be resolved by adding a few steps of (for example) Wuppertal smearing~\cite{Gusken:1989ad} as it is done in our $\Omega$-study in Chapter~\ref{cha:LatticeTechniques}.\\
However, there are already working studies on the reconstruction of the tail of the vector meson for Domain-wall fermions~\cite{RBC:2024fic}, Wilson fermions~\cite{Djukanovic:2024cmq} and also staggered fermions~\cite{Lahert:2024vvu,Lahert:2023ore}. However, in the staggered study of the MILC collaboration only a coarse lattice ($a \approx 0.15\,\mathrm{fm})$ is used with a physical extension of less than $5\,\mathrm{fm}$. This leads to the fact that only three different two-pion states are included: Goldstone taste with one or two units of momentum and the axial taste with one unit of momentum. In combination with an appropriate smearing of the vector correlator they are able to give a proper reconstruction of the long-distance regime. Other discretization schemes of fermions, on the other hand, already provide a continuum extrapolation of the long-distance reconstruction~\cite{Djukanovic:2024cmq,RBC:2024fic}.\\
In summary it can be said that the reconstruction of the tail of the vector correlator using staggered fermions is hard task and there is still a lot of work to do. For other discretizations this task is easier already applied to $g-2$. Luckily, there are also methods like data-driven approaches, that allow us to compute the long-distance contributions to $g-2$ with a high-precision~\cite{Davier:2023cyp,Boccaletti:2024guq}.