\chapter{Global analysis}
\label{cha:GlobalAnalysis}
In the following chapter, the choice of the scheme, i.e., our way to set the scale and the physical point, and our global fit procedure are described.
\section{Schemes, transformations and global fit functions}
\label{sec:scheme_fkp}
The pseudo-scalar decay constant can be decomposed in the QCD and QED contributions according to
\begin{align}
 \mathcal{F}_H^2 = f_H^2(1+\delta_H)\quad (H = \pi,K).
\end{align}
$af_H$ is a quantity that can be measured fairly easy using lattice QCD as it is shown in Section~\ref{sec:mass_fkp}. Thus, it is a useful (and common~\cite{DiCarlo:2019thl,Borsanyi:2013bia,BMW:2012hcm}) strategy to set the lattice spacing $a$ with this quantity. The problem arises that a proper physical (i.e., scheme independent) value does not exist for $f_H$ and $\delta_H$, respectively. Basically, one is free to choose how to split up the QCD and QED contributions. However, the total QCD+QED decay constant is a physical quantity that is related to the decay width by some kinematical factors. 
The lattice spacing in this study is set according to 
\begin{align}
\label{eq:set_a}
 a = \frac{a f_\pi}{f_\pi^\text{phys}} \quad \textrm{with }f_\pi^\textrm{phys} = 130.5\,\mathrm{MeV}.
\end{align}
Here, $f_\pi^\text{phys}$ is taken from the FLAG report 2024~\cite{FlavourLatticeAveragingGroupFLAG:2024oxs}. It should be mentioned that in the pure QCD limit $e \to 0$, $\mathcal{F}_H$ and $f_H$ are the same quantity and can therefore be used analogously.\\
$\mathcal{F}_\pi$ can be determined experimentally with high precision~\cite{ParticleDataGroup:2024cfk,FlavourLatticeAveragingGroupFLAG:2021npn}. In \cite{DiCarlo:2019thl} the radiative corrections are then determined using lattice simulations and the pion decay constant can thus be corrected for electromagnetic contributions. In this way, an independent value of the pure QCD pion decay constant is obtained, which is given by
\begin{align}
    f_\pi^\textrm{GRS} = 130.65(12)\,\mathrm{MeV}.
\end{align}
This quantity has a one-permille uncertainty that is highly dominated by the lattice estimation of the radiative corrections.
The ratio of the full QCD+QED pseudo-scalar decay constant can be expanded around the physical point with our ``master equation'':
\begin{align}
\begin{aligned}
\label{eq:mastereqation}
    \frac{\mathcal{F}_K}{\mathcal{F}_\pi} = A &+ B\left[\frac{M_\pi^2}{\mathcal{F}_\pi^2} - \text{const.}\right]  +C\left[\frac{M_{ds}^2 + M_{us}^2 - M_\pi^2}{\mathcal{F}_\pi^2} - \text{const.}'\right]
    \\ &+  D\left[\frac{M_{ds}^2 - M_{us}^2}{\mathcal{F}_\pi^2} - \text{const.}''\right] + Ee_v^2 + Fe_ve_s + G e_s^2.
\end{aligned}
\end{align}
Here $e_v$ denotes the charge of a valence quark and $e_s$ is the charge of a sea quark.
The constants are defined by the scheme. Here, we use the FLAG scheme~\cite{FlavourLatticeAveragingGroupFLAG:2024oxs} since the quantities are defined without any uncertainty. However, other groups~\cite{DiCarlo:2019thl} might use another scheme like the GRS, but we are able to transform our results into any scheme afterwards. The FLAG scheme is defined by
\begin{align}
    f_\pi^\textrm{FLAG} = 130.5\,\mathrm{MeV},\\
    M_\pi^\textrm{FLAG} = 135.0\,\mathrm{MeV},\\
    M_{K^+}^\textrm{FLAG} = 491.6\,\mathrm{MeV},\\
    M_{K^0}^\textrm{FLAG} = 497.6\,\mathrm{MeV}.
\end{align}
The constants $A$, $B$ and $C$ can be determined from fits to isospin-symmetric data in pure QCD
\begin{align}
\begin{aligned}
\label{eq:fkapi_iso}
    \left[\frac{f_K}{f_\pi}\right]^\textrm{isoQCD} = A + B\left[\frac{M_{ll}^2}{f_{ll}^2} - \text{const.}\right] + C\left[\frac{2M_{ls}^2- M_{ll}^2}{f_{ll}^2} - \text{const.}'\right]
\end{aligned}
\end{align}
with the constants set to the according physical values in the FLAG scheme. Since these terms are symmetric with respect to the interchange of up- and down quark we can compute the antisymmetric derivative with respect to their difference and find an expression for $D$:
\begin{align}
    D = \frac{(f_{ds} - f_{us})f_{ud}}{2(M_{ds}^2 - M_{us}^2)}.
\end{align}
This expression can be rescaled by $\left[ \frac{{M_K^0}^2 - {M_K^+}^2}{2f_\pi^2}\right]^\mathrm{FLAG}$ to get the isospin-breaking contributions to the ratio of the decay constants. The value of E can be taken from~\cite{DiCarlo:2019thl}; however, care must be taken regarding the scheme, as the GRS scheme is used there instead of the FLAG scheme. However, we can compute the transformation and the difference will be small. The $F$ term has a flavor $SU(3)$ suppression since $ \sum\limits_{s = u,d,s}  e_s = 0$ and therefore is neglected within this study. Using the notation of~\cite{Borsanyi_2021} the contribution of two sea quarks can be computed via
\begin{align}
    G = \left[\frac{\mathcal{F}_K}{\mathcal{F}\pi} - B\left[\frac{M_\pi^2}{\mathcal{F}_\pi^2} - \text{const.}\right]  -C\left[\frac{M_{ds}^2 + M_{us}^2 - M_\pi^2}{\mathcal{F}_\pi^2} \right]\right]''_{02}.
\end{align}
Here, the parameters $B$ and $C$ can be taken from the isospin-symmetric fits in Equation~\eqref{eq:fkapi_iso}. Using the parameters it is easy to transform to a given scheme by adjusting the constants to the desired quantities. For example, the PDG values can be used to obtain the physical result. The study in~\cite{DiCarlo:2019thl} from which we take the value of $E$ uses the GRS scheme~\cite{Gasser:2003hk}. In this scheme the pseudo-scalar masses are split up into their QCD and QED contributions. In these meson matching computations they neglected the sea quark effect, thus, the QED part of their calculation only consists of valence-valence effects and thus the QCD part of the scheme-defining quantities also contains sea-sea and sea-valence effects.\\
The ``master equations'' for the observable $f_\pi w_0$ and $f_\pi \sqrt{t_0}$ can be constructed in the same way. However, Equation~\eqref{eq:mastereqation} includes strong-isospin breaking effects which can be neglected here.



\section{Global fit procedure}
\label{sec:global_fit}
At this point, we have data for the isospin-symmetric and the isospin-breaking (IB) parts of $f_K/f_\pi$ at infinite volume, finite lattice spacing and quark masses scattering around the physical point. Next we want to extrapolate the data to the physical point\footnote{This is the point at which the expectation values of certain observables (e.g. pseudo-scalar masses, ...) match their physical counterparts} as well as the continuum $a = 0$. We perform three different global fits according to
\begin{align}
\label{eq:globalfit_iso}
 Y &=A +B X_l + CX_s,\\ 
\label{eq:globalfit_ib}
 D &= A_\textrm{IB} +B_\textrm{IB} X_l + C_\textrm{IB}X_s,\\
 \label{eq:globalfit_sea}
G &= A_\textrm{sea}.
\end{align}
The quantities $A$, $B$ and $C$ are computed using fits to our isospin-symmetric data, while the isospin-breaking data are fitted to $D$. The sea-quark contributions $G$ are taken from fits the QCD+QED data. $X_l$ is a quantity proportional to the light quark mass deviation from the physical point and $X_s$ is proportional to the strange quark mass deviation. In Equation~\eqref{eq:globalfit_ib} also the $X_l$ (and no $X_u$ or $X_d$) is used because this would make up second order effects in the isospin-breaking which are ignored in this study. $X_l$ and $X_s$ are used to interpolate the data to the physical point. This can be done very precisely, because the ensembles are defined in a way that they scatter closely around the physical point (see Section~\ref{sec:config})~\cite{Borsanyi_2021,Boccaletti:2024guq}. The sea quark effects are not interpolated to the physical point, since the meson masses do not scatter aroung the physical value due to large finite volume effects. However, an additional finite volume scaling is included in $A_\textrm{sea}$. The quantities $A$, $B$ and $C$ (and their isospin-breaking counterparts) are function depending on the lattice spacing that can be used to extrapolate the data to the continuum. In this study we choose this three functions to be polynomials in $a^2$ respectively $\Delta_{KS}$. The quantity $a^2$ is set using the ratio $\left(\frac{af_\pi}{f_\pi^\text{phys}}\right)^2$, where the numerator can be measured on each ensemble and the denominator can be taken from~\cite{FlavourLatticeAveragingGroupFLAG:2021npn}. Alternatively, one could use the $w_0$-scale~\cite{BMW:2012hcm,Borsanyi:2022ygn} that is used e.g. in the $g-2$ studies~\cite{Borsanyi_2021,Boccaletti:2024guq} and is further explained in Section~\ref{cha:Pheno}. $\Delta_{KS}$ is used as another parameter to quantify the discretization effects. It is much more sensitive to the lattice spacing than $a^2$, in fact it is $a^2 \alpha_s(1/a)^3$, where $\alpha_s$ is the strong coupling at the energy scale $1/a$. In Section~\ref{sec:config} it could be shown that at small lattice spacing this is close to an $a^4$ scaling. Of course, on the lattice only $a^2 \Delta_{KS}$ could be measured so the definition of $a^2$ that was mentioned before is used to remove the $a^2$ prefactor. 
\subsection{Isospin symmetric}
For the isospin symmetric fit the quantity for the light quark mass deviation is set as
\begin{align}
 X_l = \frac{M_{ll}^2}{f_{ll}^2} - \left[\frac{\left(M_{\pi^\pm}^\text{phys}\right)^2}{\left( f_{\pi^\pm}^\text{phys}\right)^2} \right]
\end{align}
and the strange quark mass deviation is given as
\begin{align}
  X_s = \frac{M_{ls}^2 - \frac{1}{2}M_{ll}^2}{f_{ll}^2} - \left[\frac{\left(M_{K^\pm}^\text{phys}\right)^2 - \frac{1}{2}\left(M_{\pi^\pm}^\text{phys}\right)^2}{\left( f_{\pi^\pm}^\text{phys}\right)^2} \right].
\end{align}
The function $A$ is a polynomial in $a^2$ and $\Delta_{KS}$ with the constraint that the sum of the maximum degrees in $a^2$ and $\Delta_{KS}$ is larger than zero and smaller than four. The function $B$ and $C$ are constructed in a similar way, but the maximum total degree is one, whereas $B$ and $C$ are also allowed to be lattice spacing independent. Since continuum extrapolation is a feature of Symanzik effective theory, which states that a lattice discretization is an asymptotic expansion in the lattice spacing around the continuum, it is possible to have a finite radius of convergence~\cite{10.1007/3-540-11192-1_11,SYMANZIK1983187,doi:10.1142/9789812777270_0004}. Therefore, it is possible that some of the ensembles with larger lattice spacings should not be included in the continuum extrapolation. In order to deal with this issue we decide to skip up to four of the largest lattice spacings. Of course, we make sure to reduce the degree of the $A/B/C$-polynomials so that over-fitting is avoided, e.g. if we remove the four largest spacing, only four spacings remain, therefore $A$ should not have more than three free parameters ($\mathrm{deg}(A) \leq 2$). We use the same argument for the $B$ and $C$ polynomial, but only the five spacing with more than two ensembles are counted.\\
There is one further exception for the polynomial allowed within the analysis:
\begin{align}
 A(a^2,\Delta_{KS}) = A_0 + A_1 a^2 + A_2a^4 + A'_1\Delta_{KS}
\end{align}
This polynomial has $a^4$ and $\Delta_{KS}$ contributions. From Figure~\ref{fig:taste_violation} it is known that both terms are very similar. In this type of fits a significant instability could be observed. So we decided to skip this function. It can be assumed that the contributions from this polynomial to the overall result are covered by the polynomials
\begin{align}
  A(a^2,\Delta_{KS}) = A_0 + A_1 a^2 + A'_1\Delta_{KS},\\
 A(a^2,\Delta_{KS}) = A_0 + A_1 a^2 + A_2a^4.
\end{align}
\begin{figure}
    \centering
    \includegraphics[width=0.49\linewidth]{Figures/7/continuum_extrapolation.pdf}
    \includegraphics[width=0.49\linewidth]{Figures/7/cumulative_density.pdf}
    \caption[Representative continuum extrapolation and histogram of the pseudo-scalar decay constant ratio.]{ \textit{Left panel:} Continuum extrapolation of the decay constant ratio of kaon and pion in the isospin-symmetric limit. The shown fit is linear in $a^2$, $\Delta_{KS}$, $X_l$ and $X_s$. The points are shifted to the physical point and are averaged if they correspond to the same lattice spacing, open points are not included in the fits. Solid lines denote the fit range, dashed lines denote regions outside of the fit range. The light blue gives the statistical deviation estimated by Jackknife resampling. \textit{Right panel:} Histogram and cumulative distribution function of the continuum extrapolated results of the pseudo-scalar decay constant ratio. The green and yellow bands give the one and two $\sigma$ neighborhoods, respectively.}
            \label{fig:cont_extrapolation}
\end{figure}
An example of one of the best isospin-symmetric fits  to $f_K/f_\pi$ is shown in the right panel of Figure~\ref{fig:cont_extrapolation}. It includes taste-breaking effects proportional to $D_{KS}$ as well as usual $a^2$ discretization effects. The shift to the physical point has no lattice spacing dependence in the example and the coarsest ensemble is excluded from the fit. However, the fit aligns well with the excluded data point.
\subsection{isospin-breaking}
For the strong-isospin-breaking contributions ($D$) only six lattice spacings, with at most two ensembles each, are available (Table~\ref{tab:configs}. The light quark mass deviation is given by
\begin{align}
 X_l &= \frac{M_{ud}^2}{f_{ud}^2} - \left[\frac{\left(M_{\pi^\pm}^\text{phys}\right)^2}{\left( f_{\pi^\pm}^\text{phys}\right)^2} \right] = \frac{M_{ll}^2}{f_{ll}^2} - \left[\frac{\left(M_{\pi^\pm}^\text{phys}\right)^2}{\left( f_{\pi^\pm}^\text{phys}\right)^2} \right] + \mathcal{O}(\delta m^2),
\end{align}
which is the same expression as for the isospin-symmetric fit up to contributions in $\delta m^2$. For the strange quark mass deviation the expression is changed a little bit to
\begin{align}
  X_s = \frac{M_{us}^2 + M_{ds}^2 - M_{ll}^2}{2f_{ll}^2} - \left[\frac{\left(M_{K^\pm}^\text{phys}\right)^2 + \left(M_{K^0}^\text{phys}\right)^2 - \left(M_{\pi^\pm}^\text{phys}\right)^2}{2\left( f_{\pi^\pm}^\text{phys}\right)^2} \right]+ \mathcal{O}(\delta m^2).
\end{align}
Since the overall number of ensembles is much smaller in this case also the number of available functions has to be reduced. The maximum total degree of of the $A_{IB}$ polynomial is reduced from three to two compared to the isospin-symmetric fits. $B_{IB}$ and $C_{IB}$ are numbers in this case with the additional constraint that at least one of them is fixed to zero. A representative fit of the isospin-breaking continuum extrapolation is shown in Figure~\ref{fig:fkp_sea_ib}.\\
Of course, all of these reductions compared to to the isospin-symmetric case will lead to a larger relative spread of mean value but since isospin-breaking effects are expected to be sub-percent, this will not affect the overall precision significantly.
\subsection{Sea quark contributions}
The quantity $G$ (defined in Equation~\eqref{eq:seasea_expansion}) that describes the sea quark contributions to $f_k/f_\pi$ is measured on the ensembles in Table ~\ref{tab:qedqcd}. Since there are only three lattice spacings available and the data have a relatively large statistical uncertainty it was made the decision to include only constant and linear functions in $a^2$ for the continuum extrapolation. Due to large finite volume effects (most of the ensembles have a box size of $\sim 3\,$fm) the pion and kaon masses are not scattered around the physical point. The extrapolation of $X_l$ and $X_s$ to zero can provide large fit instabilities. Therefore, $X_l$ and $X_s$ are not considered in this study ($B = C = 0$) and we assume that the deviation of the result by shifting the physical point is covered by the large statistical errors. Since there are two $6~\textrm{fm}$ boxes available on the coarsest lattice spacing also a finite size scaling of $\frac{1}{L^2}$ is included in the analysis. Of course, usually one does not fit an explicit dependence with two independent variables (reduced number of degrees of freedom should be at least one) but we decided to include it to somehow deal with the finite volume effects. Examples of the continuum and infinite volume extrapolations can be seen in Figure~\ref{fig:fkp_sea_ib}. The \texttt{4HEX} ensemble is expected to have smaller discretization artifacts than the other ones. Therefore, when we include it in the fit its lattice artifacts are overestimated. By repeating the analysis with the assumption that the lattice spacing of this specific ensemble is exactly zero, we also include an underestimation of the discretization effects. Thus, we include a conservative estimation of the uncertainty originating from the lattice actions.
\begin{figure}
    \centering
  \includegraphics[width = 0.55\linewidth]{Figures/7/continuum_extrapolation_ib.pdf}\\
  \includegraphics[width = 0.9\linewidth]{Figures/7/continuum_extrapolation_fkp_sea.pdf}    \caption
    [Continuum extrapolation of strong-isospin breaking and sea quark effects]{\textit{Upper panel}: Continuum extrapolation of the strong-isospin breaking contributions to $f_K/f_\pi$. The fit form is a linear function in $a^2$. All points that correspond to the same lattice spacing are averaged since no quark mass dependence is considered in this fit form. The open points denote points that are not included in the overall fit. \textit{Lower panels:} Relative sea quark contributions to $f_K/f_\pi$. The data are fitted to zero lattice spacing and infinite volume simultaneously. We present the continuum extrapolation and infinite volume and the infinite volume extrapolation at zero lattice spacing. In order to do so, the QED+QCD ensembles from Table~\ref{tab:qedqcd} are used. The red points denote the \texttt{4HEX} ensemble that is not included in the averaging procedure of measurements at same $a$ or $L$, respectively. The results corresponding to this fits are given in Table~\ref{tab:fkp}.
    }
    \label{fig:fkp_sea_ib}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{Figures/8/continuum_extrapolation_fpiw0_sea.pdf}
    \caption[Sea quark contributions to $f_\pi w_0$]{Sea quark effects contributing to $f_\pi w_0$. On the left plot the continuum extrapolation at infinite volume is given while on the right plot the infinite volume extrapolation at zero lattice spacing is shown. The strategy is the same as used for the sea quark contributions to $f_K/f_\pi$ (see Section~\ref{sec:fkapi_sib} and \ref{sec:global_fit}).  The results corresponding to this fits are given in Table~\ref{tab:fpiw0}.}
    \label{fig:fpiw0sea}
\end{figure}
\section{Further aspects of the analysis}
\label{sec:further_aspects}
In this section we want to highlight some methods and ideas that are important for this study but do not fit really into one of the other sections.
\subsection{Weighting of the systematics}
Within the previous sections a lot of different systematics were mentioned. Basically, there are three kinds of systematics:
\begin{enumerate}
 \item Systematics related to the extraction of mass and decay constants from the propagator:
 \begin{itemize}
  \item Local vs. midpoint effective mass
  \item Size of the time steps in the effective mass formula $\mathrm{d}t \in \{0.2\,\mathrm{fm},0.4\,\mathrm{fm}\} $
  \item Early vs. late fit range of the mass and decay constant plateaus
 \end{itemize}
 \item Systematics related to the NSNLO finite volume corrections
 \begin{itemize}
  \item Choice of the Hairpin parameters: $(\delta_A,\delta_V) \in \{(-1,76,0.92),(-1.87,0.00)\}$
  \item NSNLO finite volume corrections: $\mathrm{NSNLO} = \mathrm{SNLO} + \{0,2\}\times (\mathrm{NNLO} - \mathrm{NLO})$\footnote{For the isospin-breaking analysis we do not compute NNLO correction but use directly $\{0,2\}\times \mathrm{SNLO}$.}
 \end{itemize}
 \item Systematics related to the global fit procedure
 \begin{itemize}
  \item Cutting out up to four\footnote{Three for isospin-breaking} of the smallest $\beta$-values
  \item $a^2$ contributions up to third\footnote{Second for IB} order
  \item $\Delta_{KS}$ contributions up to third\footnote{Second for IB} order
  \item Nine\footnote{Three for IB} total combinations of $B$ and $C$ functions
  \end{itemize}
\end{enumerate}
As it is in done in~\cite{Borsanyi_2021} and in~\cite{Boccaletti:2024guq} the systematics related to the global fit procedure are weighted according to a modified version of the Akaike-Information-Criterion (AIC)~\cite{Akaike:1974vps,Borsanyi_2021}:
\begin{align}
 p =\mathcal{N} \exp\left(-\frac{1}{2}\left(\chi_\text{fit}^2 + 2n_\text{fit} - n_\text{data}\right) \right).
\end{align}
Here $\chi^2$ denotes the total $\chi^2$ of the global fit as described in Section~\ref{sec:global_fit}, $n_\text{fit}$ is the number of parameters that are used in this specific fit and $n_\text{data}$ is the number of data points (ensembles) incorporated in this fit. $\mathcal{N}$ is a properly chosen normalization factor. The remaining systematics are ``flat-''\-weighted, which means that for a fixed fit procedure all fit have the same weight. For example if we fit the data to 
\begin{align}
 Z = A_0 + a^2A_1 + B_0X_l + C_0X_s
\end{align}
with no $\beta$ cuts and label all systematics related the mass/decay constant extraction or FV effects by $i,j$, we get a weight for each continuum limit of
\begin{align}
 w_j^Z = \frac{p_Z^j}{\sum_Z p_Z^j\sum_i 1} \\
 p_Z^j = \exp\left(-\frac{1}{2}\left(\chi_\text{Z,j}^2 -21 \right) \right),
\end{align}
where $\chi_\text{Z,j}^2$ is the $\chi^2$ of the data to the fit model $Z$ using remaining systematics $j$.
\subsection{Blinding and cross checks}
To be unbiased during the process of analyzing the data, the kaon decay constant, $w_0$ and $\sqrt{t_0}$ are multiplied with different unknown blinding factors $\alpha^\pi_\text{blind}$, $\alpha^{w_0}_\text{blind}$ and $\alpha^{t_0}_\text{blind}$ respectively. The value of the blinding factors is not known by any of the people involved in this project. The advantage of a multiplicative blinding factor is that our relative precision is unchanged by the blinding procedure, so it can be compared to precisions of other collaboration during the process.\\
Each step of the whole analysis was performed by at least two people from the group consisting of Fabian Frech, Keith Kelley, Finn Stokes, Gen Wang and Kalman Szabo. This rigorous crosscheck procedure helps ensure the absence of inadvertent errors.