\chapter{Simulation of lattice gauge theories}
\label{cha:simulation}
The regularization with a hard cutoff is not the only advantage of discretization of spacetime. In combination with finite volume and finite time-extension, the lattice theory has a finite number of degrees of freedom in contrast to an uncountable infinite number as it is in the continuum theory. Therefore, the path-integral can be solved with a combination of analytical and numerical computations.
\section{Wick's theorem}
\label{sec:wicksthm}
One of the most important statements of the spin-statistics theorem is that it is not possible to have two fermions in the same state. Therefore fermions are described by Grassmann-valued fields $\lbrace \psi^\dagger \rbrace$ and $\lbrace \psi \rbrace$ (see Section~\ref{sec:Grassmann}), that go in the fermion bilinear (see previous chapter~\ref{cha:fermions}). The total path-integral that has to be evaluated in order to get the expectation value of an observable has the form
\begin{align}
    \langle O \rangle = \frac{1}{Z}\int [\mathcal{D}U][\mathcal{D}\chi][\mathcal{D}\bar\chi]\, O[U,\chi,\bar\chi]\exp\left(-S_g[U] - \bar \chi M \chi\right),\\
 Z = \int [\mathcal{D}U][\mathcal{D}\chi][\mathcal{D}\bar\chi]\, \exp\left(-S_g[U] - \bar \chi M \chi\right).
\end{align}
While the gauge fields $U$ consists of (at most) nine complex numbers for every link on the lattice, the Grassmann field $\bar \chi$ and $\chi$ are harder to deal with. In principle a set of $N$ Grassmann numbers can be represented using matrices of size $2^N\times 2^N$. Since $N\gg 1000$ usually it is not possible to store this representations. Luckily, the Grassmann integration can be done analytically, if a moment of a Gaussian integral has to be computed~\cite{GattringerLang}:
\begin{align}
\label{eq:wicks}
    \int [\mathcal{D}\chi][\mathcal{D}\bar\chi]\, \chi_{i_1}\chi_{i_2}...\chi_{i_k}\bar\chi_{j_1}\bar\chi_{j_2}...\bar\chi_{j_k}\exp{-\bar \chi M \chi} = \delta_{kl}\times f(M,\{i\},\{j\})\\
    f(M,\{i\},\{j\}) = (-)^k \sum_{\sigma \in S_K}\mathrm{sign}(\sigma)\left[\prod_{r = 1}^k\left(M^{-1} \right)_{i_r,j_{\sigma(r)}}\right]\times \det M
\end{align}
This is called \textit{Wick's theorem}. What finally remains is an integral with respect to the gauge links:
\begin{align}
\label{eq:PI_det}
    \langle O \rangle = \frac{1}{Z}\int [\mathcal{D}U]\, \tilde O[U]\exp\left(-S_g[U] + \log \det M[U]\right),\\
 Z = \int [\mathcal{D}U]\, \exp\left(-S_g[U] + \log \det M[U] \right),
\end{align}
where $\tilde O[U]$ is $O[U,\chi,\bar\chi]$ after expanding in powers of $\chi$ and $\bar \chi$ and applying Wick's theorem~\cite{GattringerLang,Peskin}.
\section{MCMC/Sampling of Monte-Carlo data}
\label{sec:mcmc}
The dimension $d$ of an integral like the one from Equation~\eqref{eq:PI_det} can easily exceed $10^6$. The computational cost when approximating such integrals with standard methods like numerical Riemann-Integration increases exponentially with $d$. In addition, the bulk of the phase space is suppressed by the exponential, so a huge amount of the computational time would be wasted (if it were even possible to perform this kind of integrations). Instead one can choose choose $N$ representative configurations $U_i$ from the distribution
\begin{align}
    p[U] = \frac{1}{Z}\exp\left(-S_g[U] + \log\det M[U]\right)
\end{align}
and compute the observables 
\begin{align}
    O_i = \tilde O[U_i]
\end{align}
on these representatives. The mean value of this measurements can be used as an estimator for the expectation value of the observable with an uncertainty proportional to the inverse square root of the number of configurations~\cite{GattringerLang,Caflisch_1998}:
\begin{align}
    \langle O \rangle = \frac{1}{N}\sum_{i = 1}^NO_i + \mathcal{O}\left( \frac{1}{\sqrt{N}}\right) 
\end{align}
This approach to integration is known as Monte Carlo (MC) integration. It relies on the random sampling of configurations, much like drawing numbers from a lottery, which is where the method gets its name. The hard task is of course to find a method that generates configurations according to a given and highly non-trivial distribution. For some simpler models than QCD there are methods for a direct sampling of configurations~\cite{Knechtli:2003yt}, but in QCD Markow-Chain-Monte-Carlo (MCMC) methods have to be used. These kinds of algorithms generate the configurations $U_i$ using the set $\{U_0,...,U_{i-1}\}$ in combination with some (pseudo-)random input. This of course leads to a certain amount of correlation between the configurations $U_i$ (and the respective observables $O_i$) that has to be taken into account. But this issue is explained in more detail in Section~\ref{sec:stat_err}.
The chain of the configurations $U_i$ is called a Markow-chain and has to fulfill the following four properties~\cite{Varnhorst2019lecture}:
\begin{enumerate}
    \item The probability from that the system evolves from a state $U$ to a state $U'$ ($W(U,U')$) within one update is in the interval $[0,1]$. This statement sounds trivial but is an important aspect of the mathematical formulation.
    \begin{align}
        0 \leq W(U,U') \leq 1
    \end{align}
    \item The sum of all probabilities to evolution with the same initial configuration is one:
    \begin{align}
        \int [\mathcal{D}U]\, W(U,U') = 1
    \end{align}
    \item There is a finite probability for the system to evolve from a given configuration $U$ to any other given configuration $U'$ after a certain amount of steps $k$.
    \begin{align}
        W^k(U,U') > 0
    \end{align}
    \item The transition probability $W(U,U')$ respects the desired distribution $p[U]$.
    \begin{align}
        W(U,U')p[U] = W(U',U)p[U']
    \end{align}
\end{enumerate}
The first two conditions ensure that $W(U,U')$ is a probability density function for a given $U$. The third condition is called ergodicity and makes sure, that all configurations can be sampled after a finite time and the last condition (detailed balance) takes care that the configurations $U_i$ from the Markow-Chain are samples according to the desired probability distribution $p[U]$.\\
 To start a Markow-chain, one needs an initial configurations. For this one can choose one of two methods: a cold start or a hot start. For a cold start, the initial gauge field configuration is trivial, a common choice is $U_\mu(x) = \mathbbm{1}$ for all $\mu$ and $x$. For a hot start, one can either use an already existing gauge configuration from the equilibrium of the algorithm or randomly chosen gauge links. If the initial configuration is outside of the equilibrium a thermalization is needed, which means that a certain amount of updates is performed before the measurement of the observables\footnote{that go in the final analysis of the data} can be started. To find out when the system has arrived in the equilibrium one can have look at certain observables and see at which point in the chain they start to build a plateau.\\
In the following two sections, two examples of MCMC-algorithms will be discussed: The Metropolis algorithm, which is comparatively simple to understand and implement and the Hybrid-Monte-Carlo algorithm, that is most commonly used in lattice gauge theories. Nevertheless, there are more algorithms like heat bath~\cite{PhysRevLett.43.553,PhysRevD.21.2308,FABRICIUS1984459,KENNEDY1985393}, overrelaxation~\cite{PhysRevD.23.2901,DUANE1987216,PhysRevLett.58.2394} and instanton~\cite{Durr:2012te,Eichhorn:2021ccz,Eichhorn:2022wxn} updates and modifications like Metadynamics~\cite{Eichhorn:2021ccz,Eichhorn:2022wxn,Eichhorn:2023uge,Rouenhoff:2022seh,Laio_2016,Eichhorn:2025aph} and Density of States~\cite{Fodor:2007vv,Gattringer_2019}. Comparative studies can be found in~\cite{Eichhorn:2021ccz,Eichhorn:2022wxn}.
\subsection{Metropolis algorithm}
The Metropolis algorithm was invented by Metropolis \textit{et al.} in 1953~\cite{Metropolis:1953am}. I will skip a general description of the algorithm and focus on the application to lattice gauge theories. Let $U_{i-1}$ by a given configuration with the action $S_\text{old} = S_g[U_{i-1}] - \log\det M[U]$. Then one link $U_{\nu}(y)$ of the configuration $U_{i-1}$ is chosen in addition to a random element of the gauge group $V$ that is close to the unity matrix. A proposal for a new gauge configuration $U'$ is made by 
\begin{align}
    U'_{\mu}(x) &= U_{\mu}(x)\cdot V\quad &\textrm{if $x = y$ and $\mu = \nu$,}\\
    U'_{\mu}(x) &= U_\mu(x)\quad&\textrm{else.}
\end{align} 
On the configuration $U'$ the new action $S_\text{new} = S_g[U'] - \log\det M[U']$ is measured. The configuration $U'$ is chosen to be the new configuration $U_{i}$ with the probability $\mathrm{min}\{1,e^{-(S_\text{new} - S_{\text{old}})}\}$. Otherwise, $U_i$ is set again to be $U_{i-1}$. This is called the Metropolis step. The acceptance rate of the new configuration can be tuned by choosing the width of the distribution of $V$ around the unity. Commonly this algorithm is only used if no quarks are considered, i.e., $M = 0$, because the computation of the local changes is too expensive, especially because the change of the gauge field in every step is very small, so a lot of determinants have to be computed to get new independent gauge configurations.
\subsection{Hybrid-Monte-Carlo algorithm}
To solve issue with the high cost of the simulation of quarks, Duane \textit{et al.} came up with another approach~\cite{Duane:1987de}, called Hybrid-Monte-Carlo (HMC) algorithm. If all of the links are changed by a bit in a way that does not have a huge effect on the action, a larger change of the gauge configuration can be performed before the Metropolis step and all-in-all less determinants have to be computed. The construction of the algorithm starts by bringing to mind the path-integral expectation value
\begin{align}
    \langle O \rangle = \frac{\int[\mathcal{D}U]\,O[U] \exp\left(-S[U]\right)}{\int[\mathcal{D}U]\,\exp\left(-S[U]\right)}.
\end{align}
Now the set of canonical momenta $\Pi$ of the gauge field $U$ can be included by expanding the fraction by $\int[\mathrm{D}\Pi]\, \exp\left(-\frac{1}{2}\Pi^2\right)$:
\begin{align}
    \langle O \rangle = \frac{\int[\mathcal{D}U]\,O[U] \exp\left(-S[U]\right)\int[\mathcal{D}\Pi]\, \exp\left(-\frac{1}{2}\Pi^2\right)}{\int[\mathcal{D}U]\,\exp\left(-S[U]\right)\int[\mathcal{D}\Pi]\, \exp\left(-\frac{1}{2}\Pi^2\right)}\\
    =   \frac{\int[\mathrm{D}U][\mathrm{D}\Pi]\,O[U] \exp\left(-S[U]-\frac{1}{2}\Pi^2\right)}{\int[\mathcal{D}U][\mathcal{D}\Pi]\,\exp\left(-S[U]-\frac{1}{2}\Pi^2\right)}
\end{align}
Next, the gauge fields are distributed according to a Hamiltonian $\mathcal{H} =S[U]+\frac{1}{2}\Pi^2 $ that describes a system with fields $U$ of mass $1$ in a potential $S[U]$. Therefore new gauge configurations can be can be generated by solving the Hamiltonian equations of motion
\begin{align}
\label{eq:ham_eom}
    \frac{\mathrm{d}U}{\mathrm{d}\tau} = \frac{\partial \mathcal{H}}{\partial \Pi},\quad \frac{\mathrm{d}\Pi}{\mathrm{d}\tau} = -\frac{\partial \mathcal{H}}{\partial U},
\end{align}
with an additional simulation time $\tau$, which has to be strictly distinguished from the Euclidean time $t_E$ (or $t$) that is used as a coordinate on the lattice. The procedure of an update step using the HMC algorithm is as follows:
First, for each lattice link $U_\mu(x)$ a corresponding conjugate momentum $\Pi_\mu(x)$ is chosen randomly from the distribution $\exp \left( -\frac{1}{2}\Pi^2\right)$. Then the Hamiltonian equations of motion (Equation \eqref{eq:ham_eom}) are solved using a numerical integration scheme, because this highly coupled system of partial differential equations cannot be solved analytically. To preserve detailed balance, one uses symmetric symplectic integrators, which guarantee reversibility and volume-preservation~\cite{Kennedy:2006ax}. The simplest example is the leap-frog method, where the updates are performed in the following way:
\begin{enumerate}
   \item $U_0 = \{U_\mu(x)\}$ and $\Pi_\frac{1}{2}=\{\Pi_\mu(x)\} - \frac{\epsilon}{2}\frac{\partial \mathcal{H}}{\partial U}\vert_{\tau = 0} $
   \item $U_i = U_{i-1} + \epsilon \frac{\partial H}{\partial \Pi}\vert_{\tau = i-\frac{1}{2}} $ and $\Pi_{i + \frac{1}{2}} = U_i - \epsilon\frac{\partial H}{\partial \Pi}\vert_{\tau = i-\frac{1}{2}} $ (for $i = 1,...,n-1$)
   \item $U_n =  U_{n-1} + \epsilon \Pi_{n-\frac{1}{2}}$
\end{enumerate}
This way of alternately updating the gauge fields and the canonical momenta is superior to a standard Euler-integration scheme, which has errors of $\mathcal{O}(\epsilon)$  while the leap-frog method has errors of $\mathcal{O}(\epsilon^2)$~\cite{GattringerLang}. Also, methods up to fourth order~\cite{Omelyan:2002qkh} are used to further reduce numerical deviations. Of course one can never get rid of all deviations caused by the numerical integration, so after a trajectory with $N$ steps and a step size of $\epsilon = \frac{1}{N}$ the Hamiltonian has changed by $\Delta \mathcal{H}$. Since configurations should be sampled according to $\exp\left( - \mathcal{H}\right)$, a Metropolis step can be included to set the new configuration either to $U_n$ or back to $U$~\cite{GattringerLang,DUANE1987216}.\\
There is one further difficulty that has to be discussed. When including fermions the overall action is given by
\begin{align}
    S = S_g[U] - \log \det M[U],
\end{align}
where $S_g[U]$ is a local and computational cheap function but $\det M[U]$ is not. So computing the force $\frac{\partial \mathcal{H}}{\partial U}$ in every step is still an extremely costly procedure. So one uses a trick called pseudo-fermions. As one might remember the determinant originates from a Gaussian integral of Grassmann variables
\begin{align}
    \det M = \int [\mathcal{D} \chi][\mathcal{D}\bar \chi]\, \exp\left(-\bar \chi M \chi\right).
\end{align}
A similar expression can be given by using complex numbers (if the matrix is Hermitian and positive definite):
\begin{align}
        \det M = \frac{1}{(2\pi)^n}\int [\mathcal{D} \psi][\mathcal{D}\psi^\dagger]\, \exp\left(-\psi^\dagger M^{-1} \psi\right).
\end{align}
Especially for small changes in the Matrix $M$, it could be shown that with a random vector $\psi = M \phi$ where $\phi$ is chosen from the distribution $\exp\left(-\bar \phi \phi \right)$ 
\begin{align}
    \frac{\det M [U]}{\det M[U']} = \exp\left( -\psi^\dagger(M[U] - M[U'])\psi \right)
\end{align}
is a very good estimator, whose uncertainties do not produce a significant systematical deviation in the distribution of the measurements. So before each trajectory, the $\psi$s are drawn properly and the Hamiltonian is replaced by
\begin{align}
    \mathcal{H}^\textit{pseudo} =\frac{1}{2}\Pi^2 + S_g[U] + \psi^\dagger M^{-1}\psi.
\end{align}
Of course there is still some effort necessary to compute the expression $\psi^\dagger M^{-1}\psi$ and it derivative with respect to the gauge fields in each step but this can be done using a conjugate gradient algorithm~\cite{Hestenes&Stiefel:1952} which bounds the computational costs compared to a Gaussian solver.\\
The HMC algorithm is also useful to be applied in pure gauge theory, whenever the action is modified in a way that local changes of the link variables produce global changes of the action~\cite{Eichhorn:2023uge}.
\section{Rooting}
\label{sec:rooting}
We still have to clarify the issue that the staggered Dirac operator has three unphysical zero modes. In principle they can be interpreted as flavors, but in this case all flavors need to have the same mass, and since $m_c\approx 580 \times m_u$ one might be convinced that this is no proper assumption to simulate QCD. To gain a better understanding of this let us add another set of four quarks to the theory:
\begin{align}
    S = S_g[U] + \bar \chi_1 M \chi + \bar \chi_2 M \chi.
\end{align}
After integration of the Grassmann-fields, this will change the partition function to
\begin{align}
    Z^{\textrm{8 quarks}} = \int [\mathcal{D}U]\, \exp\left(-S_g[U]\right) \times \left( \det M\right)^2.
\end{align}
It can be concluded that doubling the number of quarks of the same kind, will square the determinant in the partition function. In the same way the number of quarks can be reduced by a factor of four by using the fourth root of the determinant:
\begin{align}
    Z^{\textrm{1 quark}} = \int [\mathcal{D}U]\, \exp\left(-S_g[U]\right) \times \left( \det M\right)^{1/4}.
\end{align}
This procedure was established in~\cite{Marinari:1981qf}, highly criticized in~\cite{Creutz:2008nk}, discussed in~\cite{sharpe2006rootedstaggeredfermionsgood} and used in many cases~\cite{Borsanyi_2021,Boccaletti:2024guq}.\\
Now the problem appears that there occur powers of $M$ unequal to one in the Hamiltonian:
\begin{align}
    \mathcal{H}^\text{pseudo, rooted} =\frac{1}{2}\Pi^2 + S_g[U] + \psi^\dagger M^{-\alpha}\psi,
\end{align}
with $\alpha$ being e.g. $\frac{1}{4}$.
There are several ideas to deal with this issue. First is to use the R-algorithm, which rewrites the determinant as an effective action
\begin{align}
    \det M^\alpha = \exp\left( \alpha \tr\log M\right) = \exp\left(-S_\textit{eff}\right),
\end{align}
where the parameter $\alpha\in\mathbb{R}$ can be chosen without any constraints. The problem with this algorithm is the computation of the force especially the gradient of the effective action, that is proportional to $\tr M^{-1}\frac{\partial M}{\partial U}$ and thus requires a full inversion of the matrix $M$. This is nearly impossible due to the enourmous size of the matrices in lattice QCD. Therefore, estimations are needed that include an additional $\mathcal{O}(\epsilon)$ contributions to the Hamiltonian. Through non area-preserving integrators, this error can be reduced to $\mathcal{O}(\epsilon^2)$. In this case the detailed-balance condition is not fulfilled, so the Metropolis step is not allowed. Therefore, the measurement will have $\mathcal{O}(\epsilon)$ corrections that have to be extrapolated to 0, which is not done usually.\\
Another choice is the Polynomial HMC (PHMC)-algorithm which tries to approximate $M^\alpha$ by a polynomial $p$ in $M$:
\begin{align}
    \det M^\alpha = \int [\mathcal{D}\psi][\mathcal{D}\psi^\dagger]\,\exp\left(-\psi^\dagger\, p(M)\,\psi\right).
\end{align}
This algorithm does not have the problem of the R-algorithm but the degree $n$ of the polynomial has to be chosen at $\mathcal{O}(100-1000)$ as one needs to save at least $n/2$ vectors at once, which increases the memory cost. In addition, comparatively large rounding errors can appear.\\
The current state-of-the-art choice is the rational HMC (RHMC)-algorithm which approximates $M^\alpha$ using a rational function 
\begin{align}
    r(M) = \sum_{i = 1}^n \alpha_i(M+\beta_i)^{-1}  
\end{align}
An advantage is that the inversions for all $\beta_i$ can be done simultaneously so that no extra inversion time is required for the molecular dynamics evaluation~\cite{FROMMER_1995}. In addition, $n$ is roughly of order 20, which is much less than what is needed for polynomials~\cite{clark2006rationalhybridmontecarlo}.
\section{Smoothing techniques}
\label{sec:smoothing}
Plain gauge configurations usually provide a lot of ultra-violet noise. This means there are local fluctuations in the gauge links that can affect the condition number of the fermion matrices, the scaling of the observables at different lattice spacings, chiral symmetry- and taste-breaking effects. To reduce this noisy contributions, one can \textit{smear} the gauge fields such that the high-energy modes are reduced.\\
In addition, also the sources for fermionic observables can be smeared to change the contribution of different modes in the correlation function~\cite{Boccaletti:2024guq}.
\subsection{Gradient flow}
\label{sec:grad_flow}
Gradient flow was used first in 2009 in the context of trivializing maps~\cite{Luscher:2009eq} and applied one year later in the context of scale setting~\cite{Luscher:2010iy}. The main idea is the local minimization of the gauge action density including the flow time $t$ (not to be confused with the Euclidean time $t_E$ or the Monte-Carlo time $\tau$) as a fifth dimension and solving the partial differential equation
\begin{align}
\label{eq:grad_flow}
    \frac{\mathrm{d}V}{\mathrm{d}t} = -g_0^2 \frac{\partial S}{\partial V}V\\
    V(t = 0) = U,
\end{align}
where $\frac{\partial S}{\partial V}$ is the $\mathfrak{su}(3)$-valued derivative of the action with respect to the gauge links~\cite{Luscher:2010iy}. The energy density measured on the configurations $V$ can be computed using perturbation theory:
\begin{align}
    \langle E \rangle = \frac{3(N_c^2-1)g^2}{128\pi^2t^2}\left(1+\bar c_1 g^2 + \mathcal{O}(g^4) \right)\\
    \bar c_1 = \frac{N_c(\frac{11}{3}L + \frac{52}{9} - 3\log 3 ) - N_f (\frac{2}{3}L + \frac{4}{9}-\frac{4}{3}\log 2)}{16\pi^2}
\end{align}
$g$ is the renormalized coupling and $L = \log(8\mu^2 t) + \gamma_E$. $\mu$ is the energy scale with $\gamma_E$ being the Euler-Mascheroni constant.
Here, the energy density $E$ is defined using the clover definition of the field strength tensor~\cite{GattringerLang,Luscher:2009eq}:
\begin{align}
\label{eq:clover}
\begin{aligned}
    E = \frac{1}{V\cdot T}\sum_x F_{\mu\nu}F_{\mu\nu}\\
    F_{\mu\nu} = \frac{-\mathrm{i}}{8}(Q_{\mu\nu} - Q_{\nu\mu})\\
    Q_{\mu\nu} = P_{\mu,\nu} + P_{\mu,-\nu} + P_{-\mu,\nu} + P_{-\mu,-\nu}
\end{aligned}
\end{align}
A visualization of the clover is given in Figure~\ref{fig:Clover_term}.
 From this it could be derived that $\langle t^2E\rangle$ is independent of the gauge coupling at leading order. At higher order this invariance is broken, but the violation vanishes in the $t \to \infty$ limit. LÜscher suggested to consider the flow time $t_0$ defined as
\begin{align}
    \langle t^2 E\rangle \vert_{t = t_0} = 0.3
\end{align}
to define a quantity that proportional to $a^2$ in order to define the spacing of a lattice
 in dependence of another gauge field.\\
 This procedure is called Gradient Flow. If we choose $S[U]$ in the flow equation~\eqref{eq:grad_flow} to be the Wilson gauge action (as it was done by Lüscher), it is called Wilson flow.\\
In case of the Wilson flow (which is a very common choice) the following equation is solved:
\begin{align}
    \frac{\mathrm{d}V_\mu(n)}{\mathrm{d}t} = -g_0^2 \left(\mathcal{P}_\mathfrak{su(3) }(\mathcal{S}_\mu(n) V_\mu^\dagger (n)) \right)V_\mu(n)\\
    V_\nu(m)\vert_{t = 0} = U_\nu(m)\quad \forall \nu, m.
\end{align}
Here, $\mathcal{S}_\mu(n)$ is the sum of the staples surrounding the link $V_\mu(n)$ at the same flow time as the link $V_\mu(n)$. $\mathcal{P}_\mathfrak{su(3)}$ denotes the projection operator on the Lie algebra $\mathfrak{su(3)}$ of anti-Hermitian traceless $3\times 3$ matrices and is defined via~\cite{TimoMaster,GattringerLang}
\begin{align}
   \mathcal{P}_\mathfrak{su(3)}(A) = \frac{A-A^\dagger}{2} - \frac{1}{6}\Tr{A-A^\dagger}.
\end{align}
This formulation comes along with $\mathcal{O}(a^2)$ discretization effects that are caused by cut-off effects of the same size in the definition of the Wilson gauge action. However, in Section~\ref{sec:improvements} it was already mentioned that the gauge action can be modified by including larger loops. This was done first by Sint and Ramos, who called this modified flow the ``Zeuthen flow''~\cite{Sint:2014pip,Ramos:2015baa}.  In the studies done in this work we use three different versions of the Zeuthen flow. All of them are based on the flow equation
\begin{align}
    \frac{\mathrm{d}V_{\mu}(n)}{\mathrm{d}t} = -g_0^2\left(1+\frac{1}{12}\nabla_\mu^* \nabla_\mu \right) \frac{\partial S_\mathrm{LW}[V]}{\partial V_\mu(n)}V_\mu(n)\\
    V_\nu(m)\vert_{t = 0} = U_\nu(m)\quad \forall \nu, m
\end{align}
Here $\nabla_\mu$ and $ \nabla^*_\mu$ denote the gauge covariant forward and backward derivatives on the lattice, respectively. These operators are defined as
\begin{align}
    \nabla_\mu \phi(n) = U_\mu(n) \phi(n + \hat \mu) - \phi(n),\\
    \nabla_\mu^* \phi(n) =  \phi(n)- U_\mu^\dagger(n-\hat\mu) \phi(n -\hat \mu).
\end{align}
The three types that we call \textit{unimproved clover (uc)}, \textit{improved clover (ic)} and \textit{Symanzik (sy)} differ by their definition of the energy density. The \textit{uc} uses the clover based energy density from Equation~\eqref{eq:clover} that is also used for the Wilson flow. For the \textit{ic}, clover and plaquette field strength tensor are added in a way that removes all of the $\mathcal{O}(a^2)$ artifacts. The third type \textit{sy} is defined with the Lüscher-Weisz improved energy density~\cite{Luscher:1984xn}, thus, it also got rid of all $\mathcal{O}(a^2)$ contribution. However, it has different $\mathcal{O}(a^4)$ contributions than the improved clover type.
\begin{figure}[H]
	\centering
	\begin{tikzpicture}[scale=0.7]
		\tikzset{EmptyDotShadow/.style={draw, circle, circular drop shadow = {on layer = back}, fill = black!25, minimum size = 6pt, inner sep=0pt, fill = gray!50}}
		\foreach \x in {0,4,8}
		{
			\draw[thick] ({\x},-0.4) -- ({\x},8.4);
		}
		\foreach \y in {0,4,8}
		{
			\draw[thick] (-0.4,{\y}) -- (8.4,{\y});
		}
		\foreach \x in {0,4,8}
		{
			\foreach \y in {0,4,8}
			{
				\node[EmptyDotShadow] at (\x,\y) {};
			}
		}
		\begin{scope}[every node/.style={sloped,allow upside down}]
			\draw[very thick] (3.2,3.6) -- node {\midarrow} (0.4,3.6) -- node {\midarrow} (0.4,0.4) -- node {\midarrow} (3.6,0.4) -- node {\midarrow} (3.6,3.2);
			%\draw[very thick] (3.6,3.6) -- node {\midarrow} (0.4,3.6) -- node {\midarrow} (0.4,0.4) -- node {\midarrow} (3.6,0.4) -- node {\midarrow} (3.6,3.6);
			%\path[very thick, white] (3.2,3.6) -- (3.6,3.6);
			%
			\draw[very thick] (4.4,3.2) -- node {\midarrow} (4.4,0.4) -- node {\midarrow} (7.6,0.4) -- node {\midarrow} (7.6,3.6) -- node {\midarrow} (4.8,3.6);
			%
			\draw[very thick] (4.8,4.4) -- node {\midarrow} (7.6,4.4) -- node {\midarrow} (7.6,7.6) -- node {\midarrow} (4.4,7.6) -- node {\midarrow} (4.4,4.8);
			%
			\draw[very thick] (3.6,4.8) -- node {\midarrow} (3.6,7.6) -- node {\midarrow} (0.4,7.6) -- node {\midarrow} (0.4,4.4) -- node {\midarrow} (3.2,4.4);
		\end{scope}
        % Labels indentifying plaquettes
        \node[] at (6,6) {$P_{\mu, \nu}(n)$};
        \node[] at (6,2) {$P_{-\nu, \mu}(n)$};
        \node[] at (2,6) {$P_{\nu, -\mu}(n)$};
        \node[] at (2,2) {$P_{-\mu, -\nu}(n)$};
	\end{tikzpicture}
	\caption[Illustration of the clover term]{Illustration of the clover term used to define the field strength tensor on the lattice~\cite{FabiTimoPP}.}
	\label{fig:Clover_term}
\end{figure}

\subsection{Smearing of the gauge field}
There are different approved methods for gauge smearing like APE\footnote{Named after the APE (Array Processor Experiment) collaboration}~\cite{ALBANESE1987163}, HYP\footnote{Named after the hypercubic structure of this method}~\cite{PhysRevD.64.034504}, HEX\footnote{Short for hypercubic exponential}~\cite{Capitani_2006} and stout\footnote{``Refers to their thick-bodied nature from the large brew of paths used in their formation (coined in a Dublin public house).''~\cite{Morningstar_2004}}~\cite{Morningstar_2004}. In this thesis I will only consider the latter one, which is nothing more than a gradient flow on the Wilson action using a simple Euler integration scheme. Usually a few steps of stout smearing $n_\text{smear}$ is applied and for each step of stout-smearing the following procedure is done: For each link the weighted sum of all staples~\cite{Varnhorst2019lecture} surrounding the link is computed:
\begin{align}
\label{eq:stout_staple}
    C_\mu(x) = \sum_{\nu\neq \mu}\rho_{\mu\nu}\left(U_\nu(x)U_\mu(x+\hat \nu)U_\nu^\dagger(x+\hat \mu) - U^\dagger_\nu(x-\hat \nu)U_\mu(x-\hat\nu)U_\nu(x-\hat \nu + \hat \mu) \right).
\end{align}
In principle $\rho_{\mu\nu}$ can be chosen arbitrary but in practice it is often set constant as~\cite{Borsanyi_2021,Zhang:2009jf} 
\begin{align}
    \rho_{\mu\nu} = \rho \sim 0.1.
\end{align}
These weighted staples are multiplied on the adjoint of the gauge links
\begin{align}
    \Omega_\mu(x) = C_\mu(x) U^\dagger_\mu(x)
\end{align}
and projected on the Lie algebra of the gauge group:
\begin{align}
    Q_\mu(x) = \frac{\mathrm{i}}{2}\left( \Omega_\mu^\dagger(x) - \Omega_\mu(x)\right) - \frac{\mathrm{i}}{2N_c}\Tr\left( \Omega_\mu^\dagger(x) - \Omega_\mu(x)\right).
\end{align}
If now $U^{(n)}_\mu(x)$ is the gauge field after $n$ applications of stout smearing and $Q^{(n)}_\mu(x)$ is the algebra-element $Q_\mu(x)$ computed using $U_\mu^{(n)}(x)$, then
\begin{align}
    U^{(n+1)}_\mu(x) = \exp\left( \mathrm{i}Q^{(n)}_\mu(x) \right)U^{(n)}_\mu(x).
\end{align}

In the limit $\rho \to 0$, $n_\text{smear}\to \infty$ and $\rho n_\text{smear} = t$ the Wilson flow at flow time $t$ is restored~\cite{Nagatsuka:2023jos,Ammer:2024hqr}.\\
Usually a few steps of smearing are applied to the gauge field before computing fermion matrices to reduce taste-breaking effects (staggered fermions~\cite{Borsanyi_2021}) or chiral symmetry breaking effects (Wilson fermions~\cite{BMW:2014pzb}). Also smearing (or the gradient flow) can be applied before measuring the topological charge
\begin{align}
    Q = \frac{1}{32\pi^2}\sum_x \epsilon_{\alpha\beta\gamma\delta}F^{\alpha \beta}(x)F^{\gamma\delta}(x)
\end{align}
which reduces UV-fluctuations that dominate the result~\cite{Eichhorn:2023uge,Rouenhoff:2024geh}.
\subsection{Smearing of the fermions sources}
\label{sec:wupp_smear}
Not only the gauge fields can be smeared but also the fermion sources. According to Equation~\eqref{eq:wicks}, fermionic operators are linear combinations and product of elements of the inverse fermion matrix. Usually it is not feasible to compute $M^{-1}$ directly but what is done is to compute the matrix elements by solving the linear system
\begin{align}
    \xi = M\eta
\end{align}
for some random sources $\xi_r$ ($r = 1,...,N_\textit{source}$) that fulfill the condition~\cite{GattringerLang}
\begin{align}
    \langle \xi_r^T \xi_s\rangle = \delta_{rs}.
\end{align}
For example, wall sources, $\xi_{t_0}(x,t) = \frac{1}{N_xN_yN_z}\delta_{t,t_0}(-)^{r(x)}$ with $r(x)$ being either 0 or 1 randomly for each lattice site, but also point sources $\xi_{t_0,x_0}(x,t) \delta_{t_0,t}\delta_{x_0,x}$ can be taken into account. To change the contributions of different excitation to the measurement of the correlation function, the point source can be smeared using a three dimensional Laplace operator~\cite{Gusken:1989ad}, thus the point source is replaced by
\begin{align}
\label{eq:Wuppertal_old}
    S(x) = (D^2 + m_{sc}^2)^{-1}\delta_{x,x_0}\delta_{t,t_0},
\end{align}
where $D^2$ is the three dimensional gauge-covariant Laplace operator on the lattice and $m_{sc}$ is a parameter that can be tuned. In the first publication~\cite{Gusken:1989ad} this parameter is set such that the smearing radius is about three lattice spacings. According to the affiliation of most of the authors of the paper that introduced this smearing method, it is called Wuppertal-smearing. A definition of Wuppertal smearing that is more suitable for practical applications was given, for example, in~\cite{Boccaletti:2024guq}. For a given pseudo-fermion field $\psi$ the Wuppertal-smearing operator $\hat W$ can be applied, which is defined via:
\begin{align}
\label{eq:Wuppertal_new}
    [\hat W \psi](x) = (1-\sigma)\psi(x) + \frac{\sigma}{6}\sum_{\mu = 1,2,3}\left(U_\mu^\textrm{3d}(x)U_\mu^\textrm{3d}(x+\hat \mu)\psi(x + 2\hat \mu) + U_\mu^{\textrm{3d},\dagger}(x-\hat \mu) U_\mu^{\textrm{3d},\dagger}(x-2\hat\mu)\psi(x-2\mu)\right)
\end{align}
Here $U^{3d}$ is the gauge field $U$ after applying a few steps of three dimensional stout-smearing, i.e., $\rho_{4\mu} = \rho_{\mu 4} = 0$ in Equation~\eqref{eq:stout_staple}. $\sigma$ is a parameter that modifies the strength of the smearing procedure and similar to stout-smearing, Wuppertal-smearing can be applied a few times on the same source. In Equation~\eqref{eq:Wuppertal_new}, it can be seen that only sites that differ by an even number of spacings in all directions are included. This is done in order to not mixing contributions from different tastes.

