\chapter{Lattice techniques}
\label{cha:LatticeTechniques}
In this chapter, I provide an overview of the lattice methods and techniques employed throughout this thesis. Particular attention is given to the basic properties of our newly generated finest ensembles. In addition, the extraction of energies and decay constants from correlation functions is explained. Furthermore, I present an overview of finite-volume effects, QED corrections, and strong isospin-breaking contributions.
\section{Action and ensembles}
\label{sec:config}
Gauge configurations are generated using a Symanzik improved~\cite{Luscher:1984xn} \texttt{4stout}~\cite{Morningstar_2004} one-link staggered action. In total, there are 31 different isospin-symmetric ensembles at eight different lattice spacings. The masses are scattered around the physical point, as shown in Figure \ref{fig:landscape}. For the sea quark contributions, an additional ensemble generated with the \texttt{4HEX} staggered action~\cite{Capitani_2006} and the DBW2 gauge action~\cite{Takaishi:1996xj} is used. A big advantage of the \texttt{4HEX} action compared to the \texttt{4stout} action is that the taste breaking is much smaller, which also reduces the associated discretization effects. However, simulations with the \texttt{4HEX} action are more likely to suffer from topological freezing at fine lattice spacings~\cite{Eichhorn:2023uge}.
\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/landscape.pdf}
 \caption[Landscape plot of the ensembles]{Landscape plots of the ensembles used in this study. The pseudo-scalar masses $M_{ll}$ and $M_{ss}$ are given in units of the $w_0$-scale normalized to their physical value. Each lattice spacing is aligned with a different color and the black point marks the physical point in the isospin-symmetric limit. The errorbars on the physical point originate from uncertainties on the physical values of the pseudo-scalar masses and the $w_0$-scale. Similar figures are shown in~\cite{Boccaletti:2024guq,Risch:2024}.}
 \label{fig:landscape}
\end{figure}
\begin{table}
\centering
 \begin{tabular}{ccccccc}
	$\beta$&$a$ [fm]&$L \times T $& $am_s$&$m_s/m_l$&\#confs&SIB\\
	\hline\hline
	3.7000&0.1315&$48\times 64$&0.057291&27.899&904&$\surd$\\
	\hline
	3.7500&0.1191&$56\times 96$&0.049593&28.038&315&\\
	&&&0.049593&26.939&516&$\surd$\\
	&&&0.051617&29.183&504&$\surd$\\
	&&&0.051617&28.038&522&\\
	&&&0.055666&28.083&215&\\
	\hline
	3.7753&0.1116&$56\times 84$& 0.047615&27.843&510&$\surd$\\
	&&&0.048567&28.400&505&\\
	&&&0.046186&26.469&507&$\surd$\\
	&&&0.049520&27.852&385&\\
	\hline
	3.8400&0.0952&$64\times 96$&0.043194&28.500&510&$\surd$\\
	&&&0.043194&30.205&436&\\
	&&&0.040750&28.007&1503&\\
	&&&0.049520&27.852&500&$\surd$\\
	\hline
	3.9200&0.0787&$80\times 128$& 0.032440&27.679&506&\\
	&&&0.034240&27.502&512&$\surd$\\
	&&&0.032000&26.512&1001&\\
	&&&0.032440&27.679&327&\\
	&&&0.033286&27.738&1450&$\surd$\\
	&&&0.034240&27.502&500&\\
	\hline
	4.0126&0.0640&$96\times 144$&0.026500&27.634&446&\\
	&&&0.026500&27.124&551&$\surd$\\
	&&&0.026500&27.634&2248&\\
	&&&0.026500&27.124&1000&\\
	&&&0.027318&27.263&985&\\
	&&&0.027318&28.695&1750&\\
	\hline
	4.1479&0.0483&$128\times 192$&0.019370&27.630&1155&\\
	&&&0.019951&27.104&1605&\\
	\hline
	4.3032&0.036&$176\times 264$&0.013751&27.630&857&\\
	&&&0.014026&27.362&267&\\
	&&&0.013751&26.825&841&\\
 \end{tabular}
 \caption
 [List of ensembles]{
	Ensembles used in this study. We also give the spacing, volume and mass parameters according to the ensembles. Ensembles that are also used for strong-isospin-breaking (SIB) simulations are marked with a tick($\surd$). Further explanations on the SIB are given in Section~\ref{sec:fkapi_sib}.
 }
      \label{tab:configs}
\end{table}
In Table \ref{tab:configs} all of the isospin symmetric \texttt{4stout} ensembles are listed. The lattice spacing $a$ and the geometric extension given by $L$ and $T$ are chosen such that the physical extension of the lattice is $\sim(6^3\times 9)\,\mathrm{fm}^4$ for all of the ensembles. The charm quark mass is not given in this table, however, it is set to $11.85$ times the mass of the strange quark~\cite{Borsanyi_2021,Boccaletti:2024guq,FlavourLatticeAveragingGroupFLAG:2021npn,McNeile:2010ji}.\\
Statistical errors are estimated with a blocked Jackknife procedure~\cite{blocking} with 48 bins (see Section \ref{sec:Jackknife}). Of course, it has to be checked that the autocorrelation of the observables is small enough, such that the corresponding autocorrelation time is smaller than the block size. Otherwise the bins might be correlated and the overall error will be underestimated. One expects autocorrelation to be stronger on finer lattices~\cite{Eichhorn:2023uge}, so it is sufficient to check the autocorrelation time on the finest ensemble ($\beta = 4.3032$). The time series of the topological charge is shown in Figure~\ref{fig:top_charge_history} while the autocorrelation function and time of the pion and kaon decay constants, the squared topological charge and the Wilson-flow based $w_0$ and $t_0$ scales can be seen in Figure \ref{fig:fkapi_ACT}.
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/7/qhist.pdf}
    \caption[Time series of the topological charge for ensembles with the finest and second finest lattice spacings.]{Time series of the topological charge for ensembles with the finest and second finest lattice spacings. The latter (``Monster'') is discussed in detail in our recent publication~\cite{Boccaletti:2024guq}, while the (new) finest ensemble first appears in~\cite{FKAPI}. Therefore, some details are investigated in this work.}
    \label{fig:top_charge_history}
\end{figure}
\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/ACF.pdf}
 \caption[Autocorrelation functions of topological charge and pseudo-scalar decay constant]{Autocorrelation functions of the kaon and pion decay constant and of the squared topological charge. The autocorrelation function and time are estimated with methods explained in~\cite{Wolff:2003sm,Joswig:2022qfe}. Even though the $x$-axis is given in units of trajectories of the HMC-algorithm, the ACT is given in units of measurements which is one twentieth of ACT in units of trajectories.}
 \label{fig:fkapi_ACT}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/7/ACTscaling.pdf}
    \caption[Lattice spacing dependency of the $w_0$-scale autocorrelation time and the topological susceptibility.]{Lattice spacing dependency of the $w_0$-scale autocorrelation time and the topological susceptibility. The autocorrelation time of $w_0$ scales approximately with $1/a^2$, while the topological susceptibility approaches the limit that was estimated in \cite{Borsanyi:2016ksw}. After dividing out the taste-splitting of the taste-singlet pion mass~\cite{MILC:2010iid} the topological susceptibility is very flat in terms of the lattice spacing.}
    \label{fig:ACT_scaling}
\end{figure}
The worst autocorrelation, that is observed, is the one from the squared topological charge, the decay constants have a significant smaller autocorrelation time. Thus, it is only necessary to check whether the autocorrelation time of the squared topological charge is smaller than the block size on the finest ensemble. This condition is satisfied for all ensembles except one at $\beta = 4.3032$, which has a relatively small number of trajectories, as shown in Table~\ref{tab:configs}. Nevertheless, the autocorrelation time of the observables of interest (decay constants and gradient-flow based scales) is significant smaller. Thus, we do not have to take care of autocorrelation effects after the blocking procedure. The scaling of the autocorrelation time of $w_0$ is also considered, and it is shown that the critical scaling is proportional to $a^2$~(see Figure~\ref{fig:ACT_scaling}). The topological susceptibility approaches, in the continuum limit, the value reported in~\cite{Borsanyi:2016ksw}.
As already mentioned in Chapter~\ref{cha:fermions} different tastes come along with different hadron masses, so the pions get heavier the more the locality of the spin-taste operator is violated. In the $SU(4)$-taste-symmetry group the tastes can be labeled according $\{P,A,T,V,S\}$, which are the abbreviations for pseudo-scalar, axial-vector, tensor, vector and scalar. In combination with the pseudo-scalar spin-structure this refers to a non-locality of the spin-taste operator by 0,1,2,3 or 4 dimensions. The local operator (pseudo-scalar spin and taste structure) is the lightest and called the Goldstone taste. The mass parameters are chosen such that the masses of these mesons bracket the physical point in Figure~\ref{fig:landscape}. The taste breaking effects are quantified by
\begin{align}
 \label{eq:DKS}
 \begin{aligned}
\Delta_{KS}(\xi) = M_\pi^2(\xi) - M_\pi^2(P),\\
\delta_{KS} (\xi) = 2\left(M_\eta^2 - M_\pi^2 - \Delta_{KS}(\xi) \right).
\end{aligned}
\end{align}
Here $\xi$ denotes the different tastes ($P$,$A$,$T$,$V$ or $S$), the $\eta$-meson is the meson corresponding to the $\lambda_8$-element of the eight dimensional irrep of the $SU(3)_F$, i.e.,
\begin{align}
 \eta = \frac{u\bar u + d\bar d - 2s\bar s}{\sqrt{6}}.
\end{align}
The taste splitting $\Delta_{KS}$ is expected to scale proportional to $a^2 \alpha_s(1/a)^3$ as already discussed in~\cite{Boccaletti:2024guq}. In this study it is also shown that this behavior is somewhat similar to an $a^4$ scaling, which is also true for the (new) finest ensemble as it can be seen in Figure~\ref{fig:taste_violation}.
\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/taste_violation.pdf}
 \caption[Taste violation]{The taste-violation $\Delta_{KS}$ of two different tastes, the line shows an exact $a^4$ scaling. The Roman number labeling is done according to~\cite{Ishizuka:1993mt}.}
 \label{fig:taste_violation}
\end{figure}


\section{Gradient flow based scale setting}
In Section~\ref{sec:grad_flow} I already mentioned the scale $t_0$ that is defined using the gradient flow, more specific the Wilson flow, by figuring out the flow-time zero at which the equation $T(t_0) = 0.3$ with $T(t) = t^2 E(t)$ holds~\cite{Luscher:2009eq,Luscher:2010iy}. $t_0$ has a mass dimension of $-1$ and can be used to set the lattice spacing (the scale) according to
\begin{align}
    a = \sqrt{\frac{t_0^\text{phys}}{t_0}}
\end{align}
with $t_0^\text{phys}$ being $0.14292(104)\,\mathrm{fm}$ for $N_f = 2+1+1$ and $0.14474(57)\,\mathrm{fm}$ for $N_f = 2+1$. The values are averages from estimations of different collaborations. The averaging procedure was performed by the Flavor Lattice Averaging Group (FLAG)~\cite{FlavourLatticeAveragingGroupFLAG:2024oxs}.
In 2012 the BMW collaboration suggested to use an alternative to $t_0$ which they called the $w_0$ scale. It is defined similarly to the $t_0$ scale, but this uses the value of $T(t)$ with $t$ being the flow time, the $w_0$ scale is defined as~\cite{BMW:2012hcm}.
\begin{align}
\label{eq:w0}
  \bigg\langle t \frac{\mathrm{d}}{\mathrm{d}t}\left[t^2 E(t)\right]\vert_{t = w_0^2}\bigg\rangle = 0.3.
\end{align}
In the original publication, the following reason for this proposal was given: The old function $t^2 E(t)$ depends only on energy scales larger than $\frac{1}{\sqrt{t}}$, i.e., if effects from small values of $t$ appear that affect $T(t)$, one will find this affects all larger values of $t$. The $w_0$ scale does not have this problem because it is defined locally. By including the derivative of the defining function $T(t)$, only effects around the energy scale $\frac{1}{\sqrt{t}}$ are considered~\cite{BMW:2012hcm}. In the first publication an $N_f = 2+1$ value for the $w_0$ scale was computed:
\begin{align}
    w_0^\text{phys} = 0.1755(18)(04)\,\mathrm{fm}
\end{align}
Thus, the scale setting variable had an uncertainty of roughly one percent of its value. Current state-of-the art lattice computations on the other hand aim for even more precise results, so more accurate estimations for the scale setting variable are required.\\
Unfortunately, it is not that simple to compute a new physical value of $w_0$. Since it is a quantity with mass dimension $-1$ the continuum limit cannot be taken (in the limit $a \to 0$ the lattice quantity $w_0$ diverges). Before taking the continuum limit we need to construct a dimensionless quantity. Usually, this is done by multiplying $w_0$ by a quantity of mass dimension one, whose physical value can be measured from experiments. A good idea is to use a hadronic observable $O$ where the non-QCD contributions are under control. The quantity then $Ow_0$ has a well-defined continuum limit, and $w_0$ can be determined by dividing by the physical value of $O$, with potential additional adjustments to account for QED effects or similar corrections, if necessary.\\
In general. two methods with different advantages are used. First, there is the mass of the $\Omega^-$-baryon that is known to a very high precision from experiments ($M_{\Omega^-} = 1672.45(29)\,\mathrm{MeV}$~\cite{ParticleDataGroup:2024cfk}). The downside is that the mass extraction of the $\Omega$-baryon on the lattice is a demanding task. Second, there is the pion decay constant, $f_\pi$, that can be computed relatively easily from the pion correlator on the lattice. However, the physical value of the pion decay constant is somewhat problematic. While the pion decay width is known up to a very high precision ($2.5281(5)\times 10^{-14}\,\mathrm{MeV}$~\cite{ParticleDataGroup:2022pth}), the radiative corrections play a significant role in the overall uncertainty. Additionally, there is the still unsolved puzzle of the CKM unitarity that might be related to a wrong estimation of $\vert V_{ud}\vert$, which is also a quantity that enters our $w_0$ computation if the pion decay constant is used.

\section{Meson mass and decay constant on the lattice}
\label{sec:mass_fkp}
The pseudo-scalar mesons on the lattice are already widely discussed in Chapter~\ref{cha:fermions} and Section~\ref{sec:config}. In order not to have to deal with taste-breaking effects, we use Goldstone or pseudo-scalar taste for the meson operators.
\begin{align}
 P_{ab}(\vec x,t) = \bar \psi_a(\vec x,t) \left(\gamma_5 \otimes \xi_5\right)\psi_b(\vec x,t)
\end{align}
Here $a$ and $b$ denote the flavor content of the quarks. They are either \textit{light} ($l$) or \textit{strange} ($s$) in most of the sections. Whenever strong-isospin-breaking is mentioned $a$ and $b$ refer to \textit{up} ($u$), \textit{down} ($d$) or \textit{strange}. The correlator of the pseudo scalar is given by
\begin{align}
 G_{ab}(t) = \frac{a^3}{L^3} \sum_{\vec x} \langle P_{ab}(\vec x,t) P_{ab}^\dagger(\vec 0,0)\rangle,
\end{align}
which can also be written as
\begin{align}
\label{eq:plain_correlator}
  G_{ab}(t) = \sum_{i = 0}^\infty \left(a_i e^{-M_i t} + b_i (-)^{t}e^{-M'_it} \right).
\end{align}
The second summand is called \textit{parity partner oscillation}~\cite{Staggered_fermions} and is only non-vanishing if the two quarks have a different mass. Thus, in our case the oscillation appears for the kaon but not for the pion. 
Since in the time direction periodic boundary conditions are employed one gets propagation in backwards time direction and additionally contributions from the exponential modes that wind around the time interval $N$ times: 
\begin{align}
  G_{ab}(t) = \sum_{i = 0}^\infty \sum_{N = 0}^\infty\left[a_i \left(e^{-M_i (t+NT)} +e^{-M_i (T-t+NT)} \right) + b_i (-)^{t} \left(e^{-M'_i (t+NT)} +e^{-M'_i (T-t+NT)} \right) \right]
\end{align}
Here $T$ is the time-extend of the lattice. One usually assumes that $M_i \gg M_0$ for $i \neq 0$, therefore, all those cases are exponentially suppressed and can be ignored, since their contribution is much smaller than the statistical noise from the evaluation of the path integral. Using the closed expression for the geometric series yields
\begin{align}
 G_{ab}(t) = a_0 \frac{e^{-M_0t} + e^{-M_0(T-t)}}{1-e^{-M_0T}} + b_0 (-)^{t}\frac{e^{-M'_0t} + e^{-M'_0(T-t)}}{1-e^{-M'_0T}}\\
 \label{eq:ps_corr}
 = 2a_0 \frac{\cosh(M_0 (T/2-t))}{e^{M_0T/2}(1-e^{-M_0T})} + 2 b_0 (-)^{t}\frac{\cosh(M'_0 (T/2-t))}{e^{M'_0T/2}(1-e^{-M'_0T})}.
\end{align}
Even for the kaon no parity partner oscillation could be observed in our data. It can be assumed that $M'_0 - M_0$ is sufficiently large such that the parity partner oscillation is suppressed at the time windows considered in this study. The coefficient $a_0$ is related to the pseudo-scalar decay constant by
\begin{align}
\label{eq:relation_decay}
 f_{ab} = 2(m_a + m_b) \sqrt{\frac{a_0}{M_0^3}}\Leftrightarrow a_0 = \frac{f_{ab}^2M_0^3}{4(m_a + m_b)^2}.
\end{align}
Putting this into Equation \eqref{eq:ps_corr} yields
\begin{align}
\label{eq:correlator_f_m}
  G_{ab}(t) = \frac{f_{ab}^2M_0^3}{2(m_a + m_b)^2}\frac{\cosh(M_0 (T/2-t))}{e^{M_0T/2}(1-e^{-M_0T})}. 
\end{align}
Here $m_a$ and $m_b$ are the quark mass parameters of the flavors $a$ and $b$.
It has to be noted that the definition of the decay constant is not unique but differs according to certain conventions. The prefactor $2$ in Equation \eqref{eq:relation_decay} can be replaced by factor of $1$ or $\sqrt{2}$ depending on the physicist's choice. Nevertheless, this will not affect the ratio of kaon and pion decay constants as long as it is done consistently.\\
The pseudo-scalar correlator is measured using  our standard random source technique (a further explanation is given in Section \ref{sec:free_pion}). From the correlator we extract the effective mass $M_0$ and the effective decay constant $f_{ab}$ by using the following procedure: The mass is the curvature of the time dependence of the correlation function. It can be computed using the local mass extraction formula:
\begin{align}
\label{eq:effmasslocal}
 M^{\textrm{eff, loc}}(t)= \frac{1}{\Delta}\cosh^{-1}\frac{G(t + \Delta) + G(t-\Delta)}{2G(t)},
 \end{align}
where $\Delta$ is some quantity that is fixed in physical units. In this study we use two different values $\Delta \approx 0.2\,\textrm{fm}$ and $\Delta \approx 0.4\,\textrm{fm}$ to estimate the systematic errors from this discretization. Additionally, $\Delta$ is always chosen to be a multiple of $2a$ in order not to include artifacts from staggered oscillation.\\
Since it is also possible that our choice of the mass extraction formula includes a systematic bias, we also include another one, the midpoint mass extraction formula
\begin{align}
\label{eq:effmassmidpoint}
  M^{\textrm{eff, mid}}(t)= -\frac{1}{2\Delta}\left[\cosh^{-1}\frac{G(t + \Delta)}{G(T/2)} - \cosh^{-1}\frac{G(t - \Delta)}{G(T/2)}\right].
\end{align}

\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/eff_mass.pdf}
 \caption[Effective mass of the pion and effective decay constant of the kaon]{\textit{Left panel:} Midpoint effective mass of the pion. The shaded regions denote the plateau windows that are chosen according to Figure \ref{fig:heatmap_ps}. \textit{Right panel:} effective decay constant of the kaon, the plateau windows are shifted by $0.5\,\mathrm{fm}$ according to the study in Figure \ref{fig:heatmap_ps}. The correlation functions that belong to this data are measured on one of the ensembles with $\beta = 3.9200$.}
 \label{fig:eff_mass}
\end{figure}
An example of the effective mass according to the midpoint formula from Equation \eqref{eq:effmassmidpoint} can be seen in Figure \ref{fig:eff_mass} in the left panel. In this figure the effective mass is only shown for $1.2\,\textrm{fm} \leq t \leq 4.0\,\textrm{fm}$. In general, it can be expected that for smaller times the higher excited states from Equation \eqref{eq:plain_correlator} can be seen, therefore the effective mass does not build a plateau (within errorbars) for smaller $t$. For larger $t$ the signal-to-noise ratio of the effective mass can decrease~\cite{Detmold_2014}, Nevertheless, for the pion this problem is not as bad~\cite{signaltonoise} as e.g. for the vector correlator~\cite{Frech:2025fjp} (see Chapter \ref{cha:Exploratory}).\\
To get an expression for the effective decay constant we plug the effective mass into Equation \eqref{eq:correlator_f_m} and solve for $f$ to get
\begin{align}
 \label{eq:eff_decay_constant}
  f^{\textrm{eff}}_{ab}(t) = \sqrt{
  \frac{2G_{ab}(t)(m_a + m_b)^2}{M^{\textrm{eff}}_{ab}(t)^3}
  \frac{e^{M^{\textrm{eff}}_{ab}(t)\frac{T}{2}}(1-e^{-M^{\textrm{eff}}_{ab}(t)T})}{\cosh\!\left(M^{\textrm{eff}}_{ab}(t-\frac{T}{2})\right)}
 }.
\end{align}
An example for the effective decay constant of the kaon can be seen in the right panel of Figure \ref{fig:eff_mass}. To extract a single value for the mass and decay constant respectively, one fits a plateau to the effective (time-dependent quantity) by minimizing either the uncorrelated or the correlated $\chi^2$:
\begin{align}
\label{eq_chi_uncorr}
 \chi^2_{\text{uncorr.}}(m) &= \sum_{t = t_\text{min}}^{t_{\text{max}}} \left(\frac{m-y(t)}{\delta y(t)}\right)^2\\
 \label{eq_chi_corr}
\chi^2_{\text{corr.}}(m) &= \sum_{t^a = t_\text{min}}^{t_{\text{max}}}\sum_{t^b = t_\text{min}}^{t_{\text{max}}} (m-y(t_a)) C^{-1}\left(y(t_a),y(t_b)\right)(m-y(t_b))
\end{align}
Here $\delta y(t)$ is the standard error (Jackknife error in this study) of $y(t)$ and $C^{-1}\left(y(t_a),y(t_b)\right)$ is the inverse covariance matrix of $y$ evaluated at $y(t_a)$ and $y(t_b)$. Usually, the correlated $\chi^2$ is the more trustworthy quantity, because it respects the correlation of the effective mass/decay constant values at different time-slices. Nevertheless, there is a problem that was also found in lattice data by Ammer and Dürr in~\cite{Ammer:2019yjp} and in experimental data~\cite{DAgostini:1993arp}: The value for $m$ from the minimization of the correlated $\chi^2$ is outside the range of values it is fit to. This is most probably an effect of the covariance matrix estimation, because it is very sensitive to the noise of the Jackknife estimator for a large amount of time-slices included (see Figure \ref{fig:dagostini_condition} and Section \ref{sec:fits}). An example of this effect for our data can be seen in Figure~\ref{fig:dagostini}, further plots on this topic can be found in~\cite{Ammer:2019yjp}.
\begin{figure}
 \centering
 \includegraphics[width=0.7\linewidth]{Figures/7/dagostini.pdf}
 \caption[d'Agostini effect at the pion mass plteau fit]{Uncorrelated and correlated plateau fit for the midpoint effective mass on one of the ensembles with $\beta = 3.8400$. $\Delta$ is set to $0.4\,\mathrm{fm}$ and the fit range is chosen to be $[1.9\,\mathrm{fm},2.9\,\mathrm{fm}$]. It can be seen that the mean value of the correlated fit is above all the mean values of the data points, whereas the uncorrelated mean value is in the middle of all the data. }
 \label{fig:dagostini}
\end{figure}
\begin{figure}
 \centering
 \includegraphics[width=0.7\linewidth]{Figures/7/test_fit_correlation_condition.pdf}
 \caption[Condition numbers of the covariance matrices of the plateau fits]{The condition number of the correlation matrix of the plateau fits. One the $x$-axis we show different observables with different systematics. Red points denote fits whose mean values are outside of the spread of the mean values of the data points. Especially in the upper-right  corner one can notice a correlation between large condition numbers and the ``d'Agostini''\ effect for the finer ensembles. This plot is provided courtesy of Finn Stokes, the plot was created for an internal meeting in 2022~\cite{finnsgraphic2022}.}
 \label{fig:dagostini_condition}
\end{figure}
From this observation the decision was made to use only uncorrelated fits for the plateaus of effective masses and decay constants.
While the effective mass seems to increase for small times and to be stable for larger times, the effective decay constant seems to decrease for both small and large times. Therefore, the total result will most probably depend significantly on the choice of the range. A good strategy to pick the correct choice is needed. Here, a strategy similar the one that was established in~\cite{Boccaletti:2024guq} is used. A fit range $[t_\text{min},t_\text{max}]$ is picked and then the $Q$-value (established in~\cite{Bruno_2023}) according to the uncorrelated plateau fit on this range is computed for all of the ensembles (see Table \ref{tab:configs}). Since the $Q$-value is expected to follow a uniform distribution, we want to look for fit ranges at which this value is distributed uniformly or shifted towards higher $Q$-values. We can quantify this strategy by performing a one-sided Kolmogorov-Smirnov (KS) test~\cite{smirnov1948table,KolmogorovSmirnov1933SullaDE}, where we compare the distribution of the $Q$-values for a given range along all ensembles to the identity function. An example of this type of test applied to the $\Omega^-$ baryon mass is shown in Figure~\ref{fig:heatmap_omega}. The $P$-value of this KS-test is shown in the heat maps in Figure~\ref{fig:heatmap_ps} for the pseudo-scalar mesons with light-light and light-strange quark content in dependence of start and end point of the fit range. According to this test we picked two ranges of length $1.0\,\mathrm{fm}$ for each meson, respectively. This is done to estimate the systematic error originating from the choice of the fit range due to contamination from excited states. For the light-light mesons we pick
\begin{align}
 [t^1_\text{min},t^1_\text{max}] = [1.9\,\mathrm{fm},2.9\,\mathrm{fm}],\\
 [t^2_\text{min},t^2_\text{max}] = [2.1\,\mathrm{fm},3.1\,\mathrm{fm}]
\end{align}
and for the light-strange quark content we decided to use
\begin{align}
 [t^1_\text{min},t^1_\text{max}] = [2.4\,\mathrm{fm},3.4\,\mathrm{fm}],\\
 [t^2_\text{min},t^2_\text{max}] = [2.6\,\mathrm{fm},3.6\,\mathrm{fm}].
\end{align}
An example of the plateau fits including statistical errors is shown in Figure~\ref{fig:eff_mass}. 
\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/heatmap_ps.pdf}
 \caption[Heatmap for estimation of the fit range]{Heatmaps that describe the quality of the plateau fits in terms of the the range $[t_{\text{min}},t_{\text{max}}]$. For each interval the $Q$-value (Section \ref{sec:qval}) is measured, the color is aligned to the $P$-value of the one-sided KS-test along all ensembles. Left-hand side: Plateau fits of the pion mass. Right-hand side: Plateau fits of the kaon mass.}
 \label{fig:heatmap_ps}
\end{figure}

\section{Finite volume effects of pseudo-scalar masses and decay constants}
\label{sec:fkp_FV}
The ensembles that are given in Table~\ref{tab:configs} have a spatial extension of around $6\,\mathrm{fm}$ and a temporal extension of around $9\,\textrm{fm}$. Therefore, the measurements on our ensembles differ from the infinite volume quantity by an amount that is not known exactly, but is expected to be at the order of magnitude of our target precision (see Figure~\ref{fig:fkapi_FV}). This problem is similar to the UV cut-off through the finite lattice spacing, but while discretization artifacts have to be removed with a continuum extrapolation, finite-volume (FV) effects are easier to deal with. The main idea in this study is to estimate the difference of the finite volume and the infinite volume quantity with chiral perturbation theory ($\chi$PT) and add this to the data, so that the global fit procedure is done using infinite volume observables. This also has the advantage that fluctuations that are aligned with variations in the volumes of the different ensembles are removed completely.\\
In the literature one can find FV corrections for the pseudo-scalar masses and decay constants up to two loops or next-to-next-to leading order (NNLO): The one loop corrections (next-to leading order) can be found in a work done by Gasser and Leutwyler~\cite{Gasser:1986vb}, while Colangelo, Dürr and Hafaeli combined an existing NNLO computation of the pion scattering, that uses Lüscher's approach, to compute NNLO finite size effects of the pseudo-scalar observables~\cite{Colangelo:2005gd}. Bijnens and Rössler made a full calculation of the NNLO FV effects for the desired quantities~\cite{Bijnens:2014dea} and provided a software framework \textbf{CHIRON}~\cite{Bijnens:2014gsa} to compute the effects depending on the volume and the low-energy constants~\cite{scherer2002introductionchiralperturbationtheory} that are taken directly from our simulation combined with the $w_0$ scale from~\cite{Borsanyi:2022ygn,Borsanyi_2021,Boccaletti:2024guq}. Including only NLO effects the FV pion decay constant is smaller than the infinite volume quantity by about $0.11\%$, whereas, the difference is $0.12\%$ if also NNLO effects are included. Therefore, the decision was made to use NLO finite volume corrections and add the difference of NNLO and NLO as a systematic error source. The resulting systematic error is much smaller than our target precision of $\sim 0.1\%$.\\
Unfortunately, staggered fermions come along with several tastes of pion that are heavier than in the continuum~\cite{Staggered_fermions} as it can be seen in Equation~\eqref{eq:DKS}, which results in a reduction of finite volume effects. However, the spectrum of the $\eta$-meson is distorted in another way (also Equation~\eqref{eq:DKS}), such that some of the tastes become lighter than in the continuum. This can affect the finite volume effects in an opposite way. In the next-to leading order case these effects are already computed for the Goldstone taste by Aubin and Bernard~\cite{Aubin:2003mg,Aubin:2003uc} and for the remaining tastes by Bailey et al.~\cite{Bailey:2011bkm,Bailey:2012jy}. The formulas in the papers mentioned above depend on the taste-violation parameters $\Delta_{KS}(\xi)$ and $\delta_{KS}(\xi)$ from Equation~\eqref{eq:DKS}. The values for these quantities can be measured on our ensembles and are translated into the hairpin parameters
\begin{align}
 \delta_A = \frac{\delta_{KS}(A)}{\Delta_{KS}(A)},\\
 \delta_V = \frac{\delta_{KS}(V)}{\Delta_{KS}(V)}.
\end{align}
The Fermilab-MILC collaboration~\cite{FermilabLattice:2018zqv} provides values of
\begin{align}
 \delta_A^{\text{MILC}} = \frac{\left(r_1^2a^2\delta_A\right)^{\text{MILC}}}{\left(r_1^2a^2\Delta_A\right)^{\text{MILC}}}\bigg\vert_{a = 0.12\,\mathrm{fm}} = \frac{-0.0958(93)}{0.053983} = -1.77(17),\\
 \delta_V^{\text{MILC}} = \frac{\left(r_1^2a^2\delta_V\right)^{\text{MILC}}}{\left(r_1^2a^2\Delta_V\right)^{\text{MILC}}}\bigg\vert_{a = 0.12\,\mathrm{fm}} = \frac{0.050(24)}{0.167563} = 0.30(14).
\end{align}
From our fits (Figure \ref{fig:fkapi_FV}) we obtain the value
\begin{align}
\label{eq:HP_BMW}
 \delta_A = -1.80
\end{align}
in combination with the observation that $\delta_V$ has no significant contribution and can there be neglected ($\delta_V = 0$).
For NNLO no staggered $\chi$PT computations are available, therefore the continuum results from Bijnens~\cite{Bijnens:2014dea,Bijnens:2014gsa} are used from the difference of NNLO and NLO contributions. As mentioned before, we have a very conservative approach to estimating the systematic error for the finite volume contributions, therefore we can expect this simplification to be covered.\\

\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{Figures/7/finite_size.pdf}
 \caption[Checks on the finite size correction models]{Finite volume effects of the pion and kaon (left column and right column) mass and decay constant (top row and bottom row). The points denote the values of simulations on different volumes. They can be compared to the predictions of NLO, NNLO and NSNLO that are given by the lines.}
 \label{fig:fkapi_FV}
\end{figure}
\begin{table}
 \centering
 \begin{tabular}{c|c|c|c|c|c}
	$\beta$&$a$ [fm]&$L \times T $& $am_s$&$m_s/m_l$&\#conf\\
	\hline\hline
	3.7000&0.1315&$32\times 64$&0.057291&27.899&938\\
	  &  &$40\times 64$&0.057291&27.899&840\\
	  &  &$48\times 64$&0.057291&27.899&904\\
	  &  &$64\times 64$&0.057291&27.899&521\\
 \end{tabular}
 \caption
 [Ensembles for the finite volume study]{
	Ensembles that are used in the FV study. The third ensemble corresponds to the $(6^3\times 9)\,\mathrm{fm}^4$ extension from the ensembles in Table \ref{tab:configs} and can also be found there.)
 }
 \label{tab:configs_fs}
\end{table}
The data for pion and kaon mass and decay constant measured on the ensembles in Table~\ref{tab:configs_fs} are shown in Figure~\ref{fig:fkapi_FV} together with three different approaches to describe the finite-volume effects, i.e., NLO, NNLO and NSNLO (NNLO with NLO contributions from staggered $\chi$PT). The four panels are combined in a single total fit with four constraints being the infinite volume results. The NSNLO also uses $\delta_A$ as a fit parameter (see Equation~\eqref{eq:HP_BMW}), while $\delta_V$ is set to zero. It is evident that including staggered effects is necessary for an appropriate description of the finite volume behavior. So it is decided to use the staggered NLO to correct for finite volume effects, while the difference of usual NLO and NNLO is used to estimate the systematic error of this approach. Additionally, we use to sets of hairpin parameters
\begin{align}
 (\delta_A, \delta_V)^1 = (-1.76,0.92)\\
 (\delta_A, \delta_V)^2 = (-1.87,0.00)
\end{align}
to account for the different values of these values by BMW and MILC collaboration.
\section{Strong isospin-breaking and QED contributions}
\label{sec:fkapi_sib}
Up to this point it has been assumed that $m_u = m_d = m_l$, i.e., the up- and down quarks have the same mass, which results in an exact $SU(2)$-isospin-symmetry of QCD. Of course, this is not correct in reality, but since in this limit the simulations are much cheaper, this approximation is used for the main contribution of the results. Nevertheless, the result of an isospin-symmetric analysis differs from the physical value by the isospin-breaking contributions. For these isospin-breaking contributions new measurements are performed using the ensembles that are marked with a tick in Table~\ref{tab:configs}. All the parameters are the same as in the isospin-symmetric ensembles but the two light quark flavors are replaced by up- and down-type quarks with the masses set according to
\begin{align}
 m_u/m_d = 0.485,\\
 \frac{m_u + m_d}{2} = m_l.
\end{align}
To quantify the IB contributions we expand the QCD value of our target observable $f_K/f_\pi$ up to first order in $\delta m = m_d - m_u$:
\begin{align}
\label{eq:IB_expansion}
 \left[f_{K^\pm}/f_{\pi^\pm} \right]^{\textrm{QCD}} = \left[f_{ls}/f_{ll} \right]^{\textrm{ISO}} + \frac{\partial}{\partial\delta m}\left[f_{K^\pm}/f_{\pi^\pm} \right]^{\textrm{QCD}}\bigg\vert_{\delta m = 0}\times \delta m^{\text{phys}} + \mathcal{O}\left({\delta m^\text{phys}}^2\right)
\end{align}
Since $\delta m$ is not an observable but a scheme-dependent quantity, it is more useful to replace it by a dimensionless observable that scales proportional to the isospin-breaking, $\frac{M_{dd}^2- M_{uu}^2}{f_\pi^2} \eqqcolon \frac{\Delta M^2}{f_\pi^2} $. The difference of the squares of the pseudo-scalar masses $\Delta M^2$ is proportional to the isospin-breaking~\cite{scherer2002introductionchiralperturbationtheory}, while $f_\pi^2$ has no LO contributions in the isospin-breaking. The expansion of Equation~\eqref{eq:IB_expansion} then becomes
\begin{align}
\label{eq:IB_expansion2}
 \left[f_{K^\pm}/f_{\pi^\pm} \right]^{\textrm{QCD}} = \left[f_{ls}/f_{ll} \right]^{\textrm{ISO}} + \frac{(f_{ds}-f_{us})f_{ll}}{2\Delta M^2}\times \left[\frac{\Delta M^2}{f_\pi^2} \right]^{\text{phys}} + \mathcal{O}\left(\left[\frac{\Delta M^4}{f_\pi^4} \right]^{\text{phys}}\right)
\end{align}
with $ \left[\frac{\Delta M^2}{f_\pi^2} \right]^{\text{phys}} = 0.34850927$~\cite{FlavourLatticeAveragingGroupFLAG:2021npn}. This value does not have an error even if it is a physical value, because the scheme is set in a way that $f_{\pi}$ and $\Delta M^2$ are fixed. More on this topic can be found in Section~\ref{sec:scheme_fkp}.\\
The quantity $ \frac{(f_{ds}-f_{us})f_{ll}}{2\Delta M^2}$ is now measured on the IB ensembles, except for $f_{ll}$ which is taken from the corresponding isospin-symmetric ensemble due to the fact that is has higher statistics.\\
As it will be discussed in Section~\ref{sec:scheme_fkp} also sea-sea contributions of can be measured, according to
\begin{align}
\label{eq:seasea}
    G = \left[\frac{\mathcal{F}_K}{\mathcal{F}_\pi} - B\left[\frac{M_\pi^2}{\mathcal{F}_\pi^2} \right]  -C\left[\frac{2M_{K}^2  - M_\pi^2}{\mathcal{F}_\pi^2} \right]\right]''_{02},
\end{align}
where $B$ and $C$ describe the quark mass dependence of the target observable. Further explanations are given in Section~\ref{sec:scheme_fkp}.
Here, $\mathcal{F}_H$ denotes the QED+QCD decay constant of the hadron $H$ which is equal to $f_H$ in the limit $e \to 0$. For simplicity, here we work in the QCD isospin symmetric limit. The notation in Equation~\eqref{eq:seasea} is a short cut for the second order contribution in the charge of the sea quarks. To derive an expression for the latter equation the second order expansions in the sea quark charge have to be included:
\begin{align}
\begin{aligned}
\label{eq:seasea_expansion}
    G = \left[\frac{\mathcal{F}_K + e_s^2 \mathcal{F}_K''}{\mathcal{F}_\pi +e_s^2 \mathcal{F}_\pi''} - B\left[\frac{(M_\pi + e_s^2 M_\pi'')^2}{(\mathcal{F}_\pi+e_s^2\mathcal{F}_\pi'')^2} \right]  -C\left[\frac{2(M_{K} + e_s^2M_{K}'')^2  - (M_\pi + e_s^2 M_\pi'')^2}{(\mathcal{F}_\pi+e_s^2\mathcal{F}_\pi'')^2} \right]\right]''_{02}\\
    =\frac{ \mathcal{F}_K}{\mathcal{F}_\pi }\left[\frac{\mathcal{F_K''}}{\mathcal{F}_K} - \frac{\mathcal{F_\pi''}}{\mathcal{F}_\pi}\right] - B\left[\frac{ 2M_\pi^2}{\mathcal{F}^2_\pi }\left[\frac{M_\pi''}{M_\pi} - \frac{\mathcal{F_\pi''}}{\mathcal{F}_\pi}\right] \right]  -C\left[\frac{ 4M_K^2}{\mathcal{F}^2_\pi }\left[\frac{M_K''}{M_K} - \frac{\mathcal{F_\pi''}}{\mathcal{F}_\pi}\right]  - \frac{ 2M_\pi^2}{\mathcal{F}^2_\pi }\left[\frac{M_\pi''}{M_\pi} - \frac{\mathcal{F_\pi''}}{\mathcal{F}_\pi}\right] \right]\\
\end{aligned}
\end{align}
The corresponding quantities $\mathcal{F}_H^{('')}$ and $M_H^{('')}$ can be measured on the QCD+QED ensembles that are given in Table~\ref{tab:qedqcd}.
\begin{table}
 \centering
 \begin{tabular}{c|c|c|c|c|c}
	$\beta$&$a$ [fm]&$L \times T $& $am_s$&$m_s/m_l$&\#conf\\
	\hline\hline
	3.7000&0.1315&$24\times 48$&0.057291&27.899&726\\
	  &  &$48\times 64$&0.057291&27.899&300\\
      \hline
	  3.7753&  0.1116&$28\times 56$&0.047615&27.843&887\\
 \hline
	  3.8400&  0.0952&$32\times 64$&0.043194&28.500&1110\\
      &&&0.043194&30.205&1072\\
      &&&0.040750&28.007&1036\\
      &&&0.039139&26.893&1035\\
      \hline\hline
0.7300&0.1120&56\times 84&0.06061&33.7280&1305
      \end{tabular}
 \caption[QED+QCD ensembles.]
 {
	Ensembles that are used in the QED+QCD study. Here, dynamical QED and QCD simulations are combined, so that sea quark effects can be considered. The first eight ensembles are generated using the \texttt{4stout} action while the last one is based on the \texttt{4HEX} action.
 }
 \label{tab:qedqcd}
\end{table}
\begin{figure}
    \centering
\begin{tikzpicture}
    \begin{feynman}
        % Quark-Antiquark-Vertex für das Pion
        \vertex (a) at (0,0);
        \vertex (b) at (2,0);
        
        % Diagramm mit geschlossener Fermionenschleife
        \diagram* {
            (a) -- [fermion, half left, edge label=\(q\)] (b),
            (b) -- [fermion, half left, edge label=\(q\)] (a),
            (a) -- [boson, edge label=\(\gamma\)] (b),
        };
    \end{feynman}
\end{tikzpicture}
\hspace{3cm}
\begin{tikzpicture}
    \begin{feynman}
        % Quark-Antiquark-Vertex für das Pion
        \vertex (a) at (0,0);
        \vertex (b) at (2,0);
        \vertex (c) at (4,0);
        % Diagramm mit geschlossener Fermionenschleife
        \diagram* {
            (a) -- [fermion, half left, edge label=\(q\)] (b),
            (b) -- [fermion, half left, edge ] (a),
            (b) -- [boson, half left, edge ] (c),
            (c) -- [boson, half left, edge label=\(\gamma\)] (b),
       };
    \end{feynman}
\end{tikzpicture}
%\hspace{3cm}

\begin{tikzpicture}
    \begin{feynman}
        % Quark-Antiquark-Vertex für das Pion
        \vertex (a) at (0,0);
        \vertex (b) at (2,0);
        \vertex (c) at (4,0);
        \vertex (d) at (6,0);
        % Diagramm mit geschlossener Fermionenschleife
        \diagram* {
            (a) -- [fermion, half left, edge label=\(q\)] (b),
            (b) -- [fermion, half left, edge ] (a),
            (c) -- [fermion, half left, edge ] (d),
            (d) -- [fermion, half left, edge label=\(q\)] (c),
            (b) -- [boson, edge label=\(\gamma\)] (c),
        };
    \end{feynman}
\end{tikzpicture}
    \caption[QED sea quark diagrams.]{Sea quark diagrams of $\mathcal{O}(e^2)$ that are included in the study.}
    \label{fig:QEDgraphs}
\end{figure}
The second derivative with respect to the sea quarks charge is given by the expression
\begin{align}
    \langle O \rangle''_{02} = \bigg\langle O_0 \frac{\det  s''_2}{\det s_0}\bigg\rangle_0 - \langle O_0 \rangle_0\bigg\langle \frac{\det  s''_2}{\det s_0}\bigg\rangle_0
\end{align}
where the subscript $0$ denotes the order of QED effects of the respective quantity or the expectation value. $s$ is the fermion matrix $\prod_f (\slashed D_\text{stagg} - m_f)^{1/4}$. Within the derivation the assumption was made that 
\begin{align}
    \bigg\langle \frac{\det  s'_1}{\det s_0}\bigg\rangle_0 = 0.
\end{align}
The latter notations and calculations were published first in~\cite{Borsanyi_2021}. The diagrams entering these computations are shown in Figure~\ref{fig:QEDgraphs}. The diagram on the right has a vertex with two quarks and two photons, which is a forbidden diagram in usual QED. However, at finite lattice spacing vertices with an arbitrary number of photons can occur. This can be derived easily by expanding one of the terms appearing in the discretized fermion action:
\begin{align}
 \psi(x)D_\mu\psi(x)& =   \frac{1}{2a} \left(\psi(x) U_\mu(x)\psi(x +a \hat\mu) -\psi(x) U_\mu^\dagger(x-a\hat \mu)\psi(x -a \hat\mu) \right)\\
 &= \frac{1}{2a} \left(\psi(x) e^{\mathrm{i}aeA_\mu(x)}\psi(x + a\hat\mu) -\psi(x)  e^{-\mathrm{i}aeA_\mu(x-a\hat\mu )}\psi(x - a\hat\mu) \right)
\end{align}
When expanding the latter expression in powers of the lattice spacing it is evident that the $a^{-1}$ contributions vanish, the constant expressions reproduce the vertex from continuum QED and there are higher contributions in $(ae)^n$ that describe a vertex with $n$ photons. This also explains why there are no diagrams with more $n$-photon vertices: They are not of $\mathcal{O}(e^2)$.\\
In addition to the \texttt{4stout} ensembles, we add some of the (more precise) \texttt{4HEX} data (see Table~\ref{tab:qedqcd}) in order to achieve a reduction of statistical noise and finite volume effects of the sea-quark contributions.\\
To compute the sea-quark effects of $\mathcal{F}_H w_0$ only sea modification to Equation~\eqref{eq:seasea} are necessary. In this case these effects are described by
\begin{align}
    G = \left[\mathcal{F}_H w_0 - B\left[\frac{M_\pi^2}{\mathcal{F}_\pi^2} \right]  -C\left[\frac{2M_{K}^2  - M_\pi^2}{\mathcal{F}_\pi^2} \right]\right]''_{02}.
\end{align}

\section{\texorpdfstring{$\Omega$}{Omega} operators and correlation functions}
There are many different operators for the positive-parity, ground-state $\Omega^-$ baryon~\cite{Ishizuka:1993mt,Golterman:1984dn,Bailey:2006zn}, but in this study we will focus on three specific operators. To be able to concisely express these operators, we define the spatial covariant shift
\begin{align}
  S_i \chi_a(x) = U_i(x)\chi_a(x + \hat i)\quad(i = 1,2,3),
\end{align}
with $\chi_a$ being the quark field with flavor content $a$. Two sequential shifts are then written as $S_{\mu\nu} = S_{\mu}S_{\nu}$. If necessary, we also use an additional ``flavor''\ index $\alpha,\beta$ or $\gamma$ that is defined and further explained in~\cite{Bailey:2006zn}.
With these definitions we can define the operators:
\begin{align}
    \Omega_{\mathrm{VI}}(t) &= \sum_{x_k \% 2 = 0} \epsilon_{abc}\left(S_1\chi_aS_{12}\chi_bS_{13}\chi_c - S_2\chi_aS_{21}\chi_bS_{23}\chi_c + S_3\chi_aS_{31}\chi_bS_{32}\chi_c\right)\\
    \Omega_{\mathrm{XI}}(t) &= \sum_{x_k \% 2 = 0} \epsilon_{abc}\left(S_1\chi_aS_{2}\chi_bS_{3}\chi_c \right)\\
    % \begin{aligned}
     \Omega_{\mathrm{Ba}}(t) &= \left[2\delta_{\alpha 1}\delta_{\beta 2}\delta_{\gamma 3} - \delta_{\alpha 3}\delta_{\beta 1}\delta_{\gamma 2}-\delta_{\alpha 2}\delta_{\beta 3}\delta_{\gamma 1} + 2\delta_{\alpha 1}\delta_{\gamma 2}\delta_{\beta 3} - \delta_{\alpha 3}\delta_{\gamma 1}\delta_{\beta 2}-\delta_{\alpha 2}\delta_{\gamma 3}\delta_{\beta 1} \right]\\
     &\sum_{x_k \% 2 = 0} \epsilon_{abc}\left(S_1\chi_{a \alpha}S_{12}\chi_{b\beta}S_{13}\chi_{c \gamma} - S_2\chi_{a\alpha}S_{21}\chi_{b \beta}S_{23}\chi_{c\gamma} + S_3\chi_{a \alpha}S_{31}\chi_{b \beta}S_{32}\chi_{c \gamma}\right) \notag
    % \end{aligned}
\end{align}
These three operators differ in the tastes they couple to, which means that their ground-state masses at finite $a$ might slightly differ. Nevertheless, in the continuum the masses of these three states should be degenerate. However, we use all three of these operators to estimate the systematic uncertainty originating from the choice of the operator. To modify the contributions of the excited states, which will yield a different operator for the GEVP, we apply Wuppertal smearing (see Section~\ref{sec:wupp_smear}) to either none, one, or two of the point sources. The Wuppertal smearing parameter $\sigma$ is set to $0.5$ and the gauge fields in the Laplacian are smeared with stout smearing and a parameter of $\rho = 0.125$. The number of smearing steps for the Wuppertal and stout smearing, respectively, can be found in Table~\ref{tab:wupp_smear}. These numbers are tuned in a way that the smearing radius (in physical units) is roughly constant along all ensembles. The number of point and smeared sources can also be found in Table~\ref{tab:wupp_smear}. Since smearing reduces the signal-to-noise ratio, the number of smeared sources is greater by $\mathcal{O}(10)$. The point sources are located on a random time-slice on the corners of a cube of size $(L/2)^3$. For each of those eight points a random (canonical) unit vector in color space is chosen.
\begin{table}[h]
    \centering
    \begin{tabular}{  c || c | c || c | c | c | c | c || c  }
	$\beta$  & $N_{\rm{Wptl}}$ & $N_{\rm{3d}}$ & $t_p$ & $t_a$ & $t_b$ & range \#1 & range \#2 & \# pt, sm sources \\
	\hline
	3.7000 &  24 &  32 & 1 & 4 & 7  &  7$\dots$15 & 8$\dots$15  &   28928, 229376  \\
	\hline
	3.7500 &  30 &  40 & 1 & 4 & 7  &  8$\dots$18 & 9$\dots$18  &   66208, 530176  \\
	\hline
	3.7553 &  34 &  46 & 1 & 4 & 7  &  9$\dots$19 & 10$\dots$19 &   61024, 488192  \\
	\hline
	3.8400 &  46 &  62 & 2 & 4 & 9  & 10$\dots$20 & 11$\dots$20 &  125440, 2807552 \\
	\hline
	3.9200 &  67 &  90 & 2 & 6 & 9  & 12$\dots$25 & 13$\dots$25 &  137472, 3038720 \\
	\hline
	4.0126 & 101 & 135 & 3 & 6 & 9  & 15$\dots$30 & 16$\dots$30 &  223360, 4235520 \\
	\hline
	4.1479 & 178 & 238 & 5 & 6 & 11 & 19$\dots$40 & 21$\dots$40 & 160544, 2068736 \\
	\hline
    \end{tabular}
\caption[Parameters of the Omega analysis]{The parameters for the $\Omega$-analysis are chosen different for each $\beta$. We give the number of smearing applications for the Wuppertal smearing ($N_\mathrm{Wptl}$) and the three dimensional stout smearing ($N_\mathrm{3d}$). $t_p$ denotes the additional shift for point sources in the GEVP, $t_a$ and $t_b$ are the time slices chosen to estimate the Eigenvector in the GEVP. To fit the correlators we use two different fit ranges $\mathrm{range }\#1$ and $\mathrm{range }\#2$. The smearing radii as well as the ranges are kept constant in physical units along all ensembles. This table is previously published in~\cite{Boccaletti:2024guq}.}
\label{tab:wupp_smear}
\end{table}

\section{Extraction of the \texorpdfstring{$\Omega$}{Omega}-baryon mass}
\label{sec:omega_gevp}
To ensure that the forward and backward-propagating states have a consistent parity we apply the folding transformation
\begin{align}
    H_t \to \begin{cases}
        \frac{1}{2}\left( H_t + (-)^{t+1}H_{T-t}\right)\quad &\textrm{if $0<t<\frac{T}{2}$}\\
        H_t\quad&\textrm{$t = 0$ or $t = \frac{T}{2}$}
    \end{cases}
\end{align}
to the $\Omega^-$-propagator $H_t$.\\
As mentioned before, the three Baryon-operators $\Omega_{\mathrm{VI}}$, $\Omega_{\mathrm{XI}}$ and $\Omega_{\mathrm{Ba}}$ couple to different excitations. Even though the ground state dominates beyond a certain time separation, contributions from higher excitations will contaminate the signal and increase overall uncertainties. Therefore, we use a GEVP (Section~\ref{sec:GEVP}) to extract the different mass excitations, similar to the reconstruction of the vector correlator in Chapter~\ref{cha:Exploratory}. This method was proposed in~\cite{Aubin_2011}. In this study, we use a six dimensional GEVP that includes point sources and smeared sources due to different contributions of excited states from the smearing procedure. Additionally, certain shifts of the propagator are included; this method is called pencil-of-functions and uses shifts of the correlation function to define different operators~\cite{DeTar:2014gla} according to
\begin{align}
    \hat O_i(t) &= \hat O(t+i),\\
    \hat O^\dagger_i(t) &= \hat O^\dagger(t-i).
\end{align}
For the point sources, an additional shift $t_p$ is introduced. In this way, the amplitude of all states entering the GEVP is of the same order of magnitude. Additionally, we can reduce the contamination from excited states that play a dominant role for early time slices in combination with point sources. In total, we use different shifts of the point and four different shifts of the smeared operator. The correlation matrix from the GEVP is given by
\begin{align}
        \mathbf{H}_t =
    \left(
    \begin{array}{cc|cccc}
	H_{t+2 t_p+0}^{pp} & H_{t + 2 t_p+1}^{pp} & H_{t+t_p+0}^{ps} & H_{t+t_p+1}^{ps} & H_{t+t_p+2}^{ps} & H_{t+t_p+3}^{ps}\\
	H_{t+2 t_p+1}^{pp} & H_{t + 2 t_p+2}^{pp} & H_{t+t_p+1}^{ps} & H_{t+t_p+2}^{ps} & H_{t+t_p+3}^{ps} & H_{t+t_p+4}^{ps}\\
	\hline
	H_{t+t_p+0}^{sp}   & H_{t+ t_p+1}^{sp}    & H_{t+0}^{ss}     & H_{t+1}^{ss}     & H_{t+2}^{ss}     & H_{t+3}^{ss}    \\
	H_{t+t_p+1}^{sp}   & H_{t+ t_p+2}^{sp}    & H_{t+1}^{ss}     & H_{t+2}^{ss}     & H_{t+3}^{ss}     & H_{t+4}^{ss}    \\
	H_{t+t_p+2}^{sp}   & H_{t+ t_p+3}^{sp}    & H_{t+2}^{ss}     & H_{t+3}^{ss}     & H_{t+4}^{ss}     & H_{t+5}^{ss}    \\
	H_{t+t_p+3}^{sp}   & H_{t+ t_p+4}^{sp}    & H_{t+3}^{ss}     & H_{t+4}^{ss}     & H_{t+5}^{ss}     & H_{t+6}^{ss}
    \end{array}
    \right)\ ,
\end{align}
where the superscripts denote the respective choice of smeared (s) and unsmeared (p) source and sink. For fixed $t_a$ and $t_b$ we solve the GEVP
\begin{align}
    \mathbf{H}_{t_a}\cdot v(t_a,t_b) = \lambda(t_a,t_b)\mathbf{H}_{t_b}\cdot v(t_a,t_b).
\end{align}
The parameters $t_a$ and $t_b$ are chosen differently for the respective ensembles. The explicit values can be found in Table~\ref{tab:wupp_smear}. To reconstruct the function $\lambda(t_a-t_b) = \lambda(t_a,t_b)$ the vector-matrix-vector product
\begin{align}
    P(t;t_a,t_b) = v^\dagger(t_a,t_b)\cdot \mathbf{H}_t \cdot v(t_a,t_b) \approx \lambda(t)
\end{align}
is computed. This way we can make sure that always the correct state is chosen, especially in the long-distance regime with large statistical noise. Further explanations are given in Section~\ref{sec:GEVP}.

\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/8/omega_Beta4.1479.pdf}
    \caption[Effective mass from the $\Omega$-propagator eigenstates]{Effective mass from the $\Omega$-propagator eigenstates on the finest ensemble with $\beta = 4.1479$. The numbers are converted to physical units by choosing $a = 0.0483\,\textrm{fm}$. \textit{Upper panel:} Excitations of the $\Omega$ baryon extracted with the GEVP. The gray horizontal line corresponds to a measurement of the Belle collaboration~\cite{Yelton_2018}. \textit{Lower panel:} Ground state effective mass from the GEVP. The blue band corresponds to a single exponential fit to the ground state propagator. Figure is previously published in~\cite{Boccaletti:2024guq}.}
    \label{fig:effmass_omega}
\end{figure}
We decided to fit single-state exponential functions to the propagators instead of using plateau fits to the effective mass. However, for better visualization, we still show the effective masses for one of the finest ensembles in Figure~\ref{fig:effmass_omega} along with the result of our ground state fit and an experimental value~\cite{Yelton_2018} for the first excitation that is compatible with our data. We fit the data to the model
\begin{align}
    Y(t) = A \exp(-Mt)
\end{align}
by minimizing the correlated $\chi^2$. The optimization with respect to the coefficient $A$ can be done analytically, so we end up with a one-parameter fit
\begin{align}
    \Tilde{Y}(t) = \frac{\sum_{t_0,t_1} \lambda(t_0) C^{-1}(t_0,t_1) \exp(-Mt_1)}{\sum_{t_0,t_1} \exp(-Mt_0) C^{-1}(t_0,t_1) \exp(-Mt_1)} \exp(-Mt),
\end{align}
where $C^{-1}(t_0,t_1)$ is the inverse covariance matrix of $y$, evaluated at $y(t_0)$ and $y(t_1)$.
Usually, we use 48 Jackknife samples to estimate the uncertainty. However, for the stability of the covariance matrix, it is useful to use a larger amount of samples. Here, we use 200 Jackknife samples, which still ensures that the blocks are large enough to handle autocorrelation effects. The autocorrelation function of the $\Omega$ propagator can be found in Figure~\ref{fig:AC_omega}.
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/8/ACF_MO1.pdf}
    \caption[Autocorrelation functions of the $\beta = 4.1479$ ensemble]{Autocorrelation functions on the $\beta = 4.1479$ ensemble. \textit{Left panel}: Autocorrelation functions of the squared topological charge and the energy density, both at flow time $w_0^2$. \textit{Right panel}: Autocorrelation function of the $\Omega_\mathrm{VI}$-operator with two point sources and two smeared sources, respectively. Figure is previously published in~\cite{Boccaletti:2024guq}.}
    \label{fig:AC_omega}
\end{figure}
The covariance matrix is inverted using an SVD-cut of $10^{-8}$, as explained in Section~\ref{sec:chi2}, to prevent the overall result from being dominated by statistical noise~\cite{Bouchard_2014}.\\
Usually, it is a non-trivial task to find the optimal region in which the exponential function is fitted to the correlator. Here, we decide to use a method very similar to the one used in Section~\ref{sec:mass_fkp}. For each ensemble, we compute the (correlated) Q-value (Section~\ref{sec:qval})~\cite{Bruno_2023} for a given range $[t_\mathrm{min},t_\mathrm{max}]$. Then, we perform a one-sided KS-test across the ensembles, where we allow the Q-values to be tilted towards 1. The P-values of this KS-tests are shown in the right panel of Figure~\ref{fig:heatmap_omega}. Based on this figure, the decision was made to use the fit ranges $[0.9~\mathrm{fm},2.0~\mathrm{fm}]$ and $[1.0~\mathrm{fm},2.0~\mathrm{fm}]$. We use two different ranges to account for the systematic uncertainty associated with different excitations. The KS-tests for both of the chosen ranges are shown in the left panel of Figure~\ref{fig:heatmap_omega}.
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/8/heatmap_omega.pdf}
    \caption[KS-test of the Q-values for different fit ranges]{\textit{Left panel:} KS-test of the Q-values along all ensembles for the ranges $[0.9~\mathrm{fm},2.0~\mathrm{fm}]$ and $[1.0~\mathrm{fm},2.0~\mathrm{fm}]$. \textit{Right panel:} Heatmap of the P-values of all KS-tests. The red crosses mark the fit ranges that are used within the final analysis. Figure is previously published in~\cite{Boccaletti:2024guq}.}
    \label{fig:heatmap_omega}
\end{figure}
All in all we are able to measure the mass of the $\Omega$-baryon on our lattices with permille precision as it can be seen in Figure~\ref{fig:global_omega}.
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/8/global_omega.pdf}
    \caption[Ground state $\Omega$ baryon mass for different ensembles]{Ground state $\Omega$ baryon mass for different ensembles. For each $\beta$ different colors denote different ensembles. The six data points per ensemble correspond to two different fit ranges times three different operators. Figure is previously published in~\cite{Boccaletti:2024guq}.}
    \label{fig:global_omega}
\end{figure}


