# Novel approach for computing gradients of physical observables

Simone Bacchio  Computation-based Science and Technology Research Center, The
Cyprus Institute, Nicosia, Cyprus

(August 28, 2024)

###### Abstract

We show that an infinitesimal step of gradient flow can be used for defining a
novel approach for computing gradients of physical observables with respect to
action parameters. Compared to the commonly used perturbative expansion, this
approach does not require calculating any disconnected contribution or vacuum
expectation value and can provide results up to three orders of magnitudes
more precise. On the other hand, it requires a non-trivial condition to be
satisfied by the flow action, the calculation of its force and its Laplacian,
and the force of the observable, whose gradient needs to be measured. As a
proof of concept, we measure gradients in $\beta$ of Wilson loops in a four-
dimensional SU(3) Yang-Mills theory simulated on a $16^{4}$ lattice using the
Wilson action.

††preprint: DESY 22-

##  I Introduction

In a lattice gauge theory with an action $S_{\theta}\equiv S(U,\theta)$, the
expectation value of an observable $\mathcal{O}(U)$ is given by the path
integral

$\displaystyle\langle\mathcal{O}\rangle_{\theta}=\frac{1}{\mathcal{Z}_{\theta}}\int\textrm{D}[U]\,\mathcal{O}(U)\,\exp(-S(U,\theta))\,,$
(1)

where $\mathcal{Z}_{\theta}\equiv\int\textrm{D}[U]\,\exp(-S(U,\theta))$ is the
partition function and $\theta$ is a parameter of the action. The focus of
this work is the calculation of the gradient of the expectation value with
respect to the parameter $\theta$, which might result in a challenging task
since the so-called sea effects have to be measured. Indeed, it is standard to
perturbatively expand the expectation value with respect to an infinitesimal
change of the parameter, $d\theta$, obtaining

$\displaystyle\langle\mathcal{O}_{\theta+d\theta}\rangle_{\theta+d\theta}=\frac{\langle
d\theta\frac{\partial\mathcal{O}_{\theta}}{\partial\theta}+\mathcal{O}_{\theta}\exp(-d\theta\frac{\partial
S_{\theta}}{\partial\theta})\rangle_{\theta}}{\langle\exp(-d\theta\frac{\partial
S_{\theta}}{\partial\theta})\rangle_{\theta}}+O(d\theta^{2})$

and, thus, the gradient with respect to $\theta$ is given by

$\displaystyle\frac{d\langle\mathcal{O}_{\theta}\rangle_{\theta}}{d\theta}$
$\displaystyle=\lim_{d\theta\rightarrow
0}\frac{\langle\mathcal{O}_{\theta+d\theta}\rangle_{\theta+d\theta}-\langle\mathcal{O}_{\theta}\rangle_{\theta}}{d\theta}$
(2)
$\displaystyle=\langle\frac{\partial\mathcal{O}_{\theta}}{\partial\theta}-\mathcal{O}_{\theta}\frac{\partial
S_{\theta}}{\partial\theta}\rangle_{\theta}+\langle\mathcal{O}_{\theta}\rangle_{\theta}\langle\frac{\partial
S_{\theta}}{\partial\theta}\rangle_{\theta}\,,$

where, in the above equations, for generality’s sake, we have introduced an
optional explicit dependence on $\theta$ also in the observable
$\mathcal{O}_{\theta}\equiv\mathcal{O}(U,\theta)$. Since on both sides of Eq.
(2) all terms depend on $\theta$, we rewrite the equation with the shorthand
notation

$\displaystyle\frac{d\langle\mathcal{O}\rangle}{d\theta}$
$\displaystyle=\langle\frac{\partial\mathcal{O}}{\partial\theta}-\mathcal{O}\frac{\partial
S}{\partial\theta}\rangle+\langle\mathcal{O}\rangle\langle\frac{\partial
S}{\partial\theta}\rangle\,$ (3)

and we refer to

$\displaystyle\langle\frac{\partial\mathcal{O}}{\partial\theta}\rangle\,:$ as
connected contribution, $\displaystyle\langle\mathcal{O}\frac{\partial
S}{\partial\theta}\rangle\,:$ as disconnected contribution, and
$\displaystyle\langle\mathcal{O}\rangle\langle\frac{\partial
S}{\partial\theta}\rangle\,:$ as vacuum expectation value.

Disconnected contributions are well-known to be noisy. Furthermore, when the
vacuum expectation value is non-zero, the latter and the disconnected
contributions are usually large values, whose difference needs to be taken
accurately for reliably measuring the gradient.

In this work, we present a novel and alternative approach, based on gradient-
flow techniques [1], which is free from any of the aforementioned problems.
Namely, if a flow action that flows along the parameter $\theta$ is found,
then gradients of observables can be computed directly using the ODE of the
operator. Our numerical results show that sea effects can be computed up to
three orders of magnitude more precisely than the perturbative expansion in
Eq. (3), when an exact flow action is found. Therefore, this approach can
possibly be used to improve the accuracy of gradients in, for example, among
others:
i.) applications of the Feynman-Hellmann theorem [2, 3, 4, 5, 6, 7],
ii.) leading isospin breaking corrections [8, 9, 10, 11, 12],
iii.) leading QED corrections [13, 14, 12],
iv.) fine-tuning of simulation parameters [15],
v.) leading contribution of the QCD $\Theta$-term to the neutron electric
dipole moment [16, 17, 18].

In the following, we present the novel approach in the form of a theorem,
followed by its proof and a numerical case study.

##  II Main result

### Notation.

In the following,
• $\tilde{S}$ is referred to as flow action and is to-be-determined;
• $\mathcal{L}_{0}=-\sum_{x,\mu,a}\partial_{x,\mu}^{a}\partial_{x,\mu}^{a}$
denotes the Laplacian;
• $(A,B)=\sum_{x,\mu,a}A^{a}_{\mu}(x)B^{a}_{\mu}(x)$ denotes the scalar
product over algebra-valued fields; and
• $\partial$ is the force of a scalar function, defined as

$\displaystyle\partial^{a}_{x,\mu}f(U)=\left.\frac{d}{d\tau}f(U_{\tau})\right|_{\tau=0}$
(4)

with

$\displaystyle U_{\tau}(y,\nu)=\begin{cases}e^{\tau T^{a}}U(x,\mu),&\text{for
}(x,\mu)=(y,\nu)\,,\\\ U(y,\nu)&\text{for }(x,\mu)\neq(y,\nu)\,.\end{cases}$
(5)

### Theorem, gradient-flow approach.

The gradient of an expectation value with respect to a parameter $\theta$ of
the action is given by

$\displaystyle\frac{d\langle\mathcal{O}\rangle}{d\theta}$
$\displaystyle=\langle\frac{\partial\mathcal{O}}{\partial\theta}+\left(\partial\mathcal{O},\partial\tilde{S}\right)-\mathcal{O}\mathcal{C}\rangle+\langle\mathcal{O}\rangle\langle\mathcal{C}\rangle\,,$
(6)

where

$\displaystyle\mathcal{C}=\mathcal{L}_{0}\tilde{S}+\left(\partial
S,\partial\tilde{S}\right)+\frac{\partial S}{\partial\theta}\,.$ (7)

### Corollary, ideal gradient-flow approach.

If the flow action $\tilde{S}$ is such that $\mathcal{C}=$ constant, then Eq.
(6) simplifies to

$\displaystyle\frac{d\langle\mathcal{O}\rangle}{d\theta}$
$\displaystyle=\langle\frac{\partial\mathcal{O}}{\partial\theta}+\left(\partial\mathcal{O},\partial\tilde{S}\right)\rangle\,.$
(8)

We, therefore, refer as ideal flow action to an $\tilde{S}$ such that

$\displaystyle\mathcal{L}_{0}\tilde{S}+\left(\partial
S,\partial\tilde{S}\right)+\frac{\partial S}{\partial\theta}={\rm
constant}\,,$ (9)

### Corollary, standard perturbative approach.

If $\tilde{S}\equiv$ constant, then Eq. (6) simplifies to Eq. (3), since
$\mathcal{L}_{0}\tilde{S}=\partial\tilde{S}=0$ and $\mathcal{C}=\partial
S_{\theta}/\partial\theta$.

### Summary.

The corollaries are two special applications of the theorem and our numerical
results show that gradients computed using the ideal gradient-flow approach
are significantly more precise than those computed using the standard
perturbative approach. If Eq. (9) is not satisfied exactly, then the theorem
can be used to obtain correct and possibly improved results.

##  III Proof

The proof is based on gradient-flow techniques that have been introduced in
Ref. [1], which we refer to for further details. Under a change of variables
$U=\mathcal{F}(V)$, the expectation value of an observable transforms as

$\displaystyle\langle\mathcal{O}\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}}\int\textrm{D}[U]\,\mathcal{O}(U)\,\exp(-S(U))$
(10)
$\displaystyle=\frac{1}{\mathcal{Z}_{\mathcal{F}}}\int\textrm{D}[V]\,\mathcal{O}(\mathcal{F}(V))\,\exp(-S_{\mathcal{F}}(V))\,$

where
$\mathcal{Z}_{\mathcal{F}}\equiv\int\textrm{D}[V]\exp(-S_{\mathcal{F}}(V))$
and

$\displaystyle S_{\mathcal{F}}(V)$
$\displaystyle=S(\mathcal{F}(V))-\ln\det\mathcal{F}_{*}(V)$ (11)

with $\mathcal{F}_{*}$ being the Jacobian of the transformation. In the
following, we consider an infinitesimal transformation of the field employing
an Euler integration step and the force of an action $\tilde{S}$

$\begin{aligned} \mathcal{F}(V)&=\exp(-\epsilon\,\partial\tilde{S}(V))V\\\
&=V-\epsilon\,\partial\tilde{S}(V)V+O(\epsilon^{2})\end{aligned}\quad\text{with}\quad
0\\!<\\!\epsilon\\!\ll\\!1.$ (12)

Two properties of such flow are [1]

$\ln\det\mathcal{F}_{*}(V)=\epsilon\,\mathcal{L}_{0}\tilde{S}(V)+O(\epsilon^{2})$
(13)

and

$\mathcal{O}(\mathcal{F}(V))=\mathcal{O}(V)-\epsilon\,\left(\partial\mathcal{O}(V),\partial\tilde{S}(V)\right)+O(\epsilon^{2})\,.$
(14)

By using Eqs. (13) and (14) in Eq. (11), we obtain

$\displaystyle S_{\mathcal{F}}=S-\epsilon\left[\left(\partial
S,\partial\tilde{S}\right)+\mathcal{L}_{0}\tilde{S}\right]+O(\epsilon^{2}),$
(15)

where all terms now depend on the field $V$. The intent is then to interpret
the difference $S_{\mathcal{F}}-S$ as an infinitesimal change of action along
the parameter $\theta$, requiring

$\displaystyle\frac{\partial S}{\partial\theta}=\lim_{\epsilon\rightarrow
0}\frac{S_{\mathcal{F}}-S}{\epsilon}=-\left(\partial
S,\partial\tilde{S}\right)-\mathcal{L}_{0}\tilde{S}\,.$ (16)

We then define $\mathcal{C}$, in Eq. (7), as the difference between the left-
and right-hand sides. By using $\mathcal{C}$ in Eq. (15), we can remove any
dependence on $\tilde{S}$ obtaining

$\displaystyle S_{\mathcal{F}}=S+\epsilon\left(\frac{\partial
S}{\partial\theta}-\mathcal{C}\right)+O(\epsilon^{2})\,.$ (17)

In the following, we track explicitly the dependence on the parameter
$\theta$, and from Eq. (17) we obtain

$\displaystyle
S_{\mathcal{F}}=S_{\theta+\epsilon}-\epsilon\,\mathcal{C_{\theta}}+O(\epsilon^{2})\,.$
(18)

By using the latter in the definition of $\mathcal{Z}_{\mathcal{F}}$, we
obtain

$\displaystyle\mathcal{Z}_{\mathcal{F}}=\mathcal{Z}_{\theta+\epsilon}-\epsilon\,\mathcal{Z}_{\theta}\langle\mathcal{C_{\theta}}\rangle_{\theta}+O(\epsilon^{2})\,.$
(19)

Finally, by using the above in Eq. (10), we obtain

$\displaystyle\langle\mathcal{O}_{\theta}\rangle_{\theta}=\langle\mathcal{O}_{\theta}\rangle_{\theta+\epsilon}$
$\displaystyle+\epsilon\,\langle\mathcal{O}_{\theta}\mathcal{C}_{\theta}\rangle_{\theta}-\epsilon\,\langle\mathcal{O}_{\theta}\rangle_{\theta}\langle\mathcal{C}_{\theta}\rangle_{\theta}$
(20)
$\displaystyle-\epsilon\,\langle\left(\partial\mathcal{O}_{\theta},\partial\tilde{S}_{\theta}\right)\rangle_{\theta}+O(\epsilon^{2})$

and the theorem is proved since

$\displaystyle\frac{d\langle\mathcal{O}_{\theta}\rangle_{\theta}}{d\theta}$
$\displaystyle=\lim_{\epsilon\rightarrow
0}\frac{\langle\mathcal{O}_{\theta+\epsilon}\rangle_{\theta+\epsilon}-\langle\mathcal{O}_{\theta}\rangle_{\theta}}{\epsilon}$
(21)
$\displaystyle=\langle\frac{\partial\mathcal{O}_{\theta}}{\partial\theta}+\left(\partial\mathcal{O}_{\theta},\partial\tilde{S}_{\theta}\right)-\mathcal{O}_{\theta}\mathcal{C}_{\theta}\rangle_{\theta}+\langle\mathcal{O}_{\theta}\rangle_{\theta}\langle\mathcal{C}_{\theta}\rangle_{\theta}\,.$

##  IV Case study and numerical results

Figure 1: Representation of the loops entering the $\beta$-expansion up to the
NLO. The plaquette, shown by the green loop, enters the LO expansion, while
the loops shown in blue enter the NLO. The enumeration of the loops follows
Ref. [1].

### Wilson action.

To demonstrate the effectiveness of Eq. (6) compared to Eq. (3), we consider a
four-dimensional SU(3) Yang-Mills theory defined on a lattice using the
standard Wilson action

$S_{W}(\beta,U)=-\frac{\beta}{6}\mathcal{W}_{0}(U),$ (22)

where $\mathcal{W}_{0}$ denotes the sum of plaquettes. We aim at computing the
slope in $\beta$ of various observables and, therefore, the approach outlined
here requires a flow action $\tilde{S}$ that satisfies Eq. (9) in the
following way,

$\mathcal{L}_{0}\tilde{S}-\frac{\beta}{6}\left(\partial\mathcal{W}_{0},\partial\tilde{S}\right)-\frac{1}{6}\mathcal{W}_{0}=\text{
constant}\,.$ (23)

Inspired by Lüscher’s $t$-expansion [1], it is easy to show that an analytical
solution to Eq. (23) is

$\displaystyle\tilde{S}$
$\displaystyle=\frac{1}{6}\sum_{k=0}^{\infty}\left(\frac{\beta}{6}\right)^{k}\tilde{S}^{(k)}\qquad\leavevmode\nobreak\
\leavevmode\nobreak\ \,\text{with}$ (24) $\displaystyle\tilde{S}^{(0)}$
$\displaystyle=\mathcal{L}_{0}^{-1}\mathcal{W}_{0}=\frac{3}{16}\mathcal{W}_{0}\quad\quad\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm and}$
$\displaystyle\tilde{S}^{(k)}$
$\displaystyle=\mathcal{L}_{0}^{-1}\left(\partial\mathcal{W}_{0},\partial\tilde{S}^{(k-1)}\right)\quad{\rm
for\leavevmode\nobreak\ }k>0\,,$

which we refer to as $\beta$-expansion and it is equivalent to Lüscher’s
$t$-expansion evaluated at $t=\beta$ divided by $\beta$. We refer to
$\tilde{S}^{(0)}$ as the leading order (LO) and to $\tilde{S}^{(1)}$ as the
next-to-leading order (NLO). The calculation of the NLO can be found in Ref.
[1], where the term
$\left(\partial\mathcal{W}_{0},\partial\mathcal{W}_{0}\right)$ is computed. It
results in a linear combination of the Wilson loops depicted in Fig. 1, which
are defined as

$\displaystyle\mathcal{W}_{i}$ $\displaystyle=\sum_{C\in\Gamma_{i}}{\rm
tr}\\{U(C)\\}$ $\displaystyle{\rm for\leavevmode\nobreak\ }i=0,1,2,5$ (25)
$\displaystyle\mathcal{W}_{i}$
$\displaystyle=\\!\\!\\!\sum_{C,C^{\prime}\in\Gamma_{i}}\\!\\!\\!{\rm
tr}\\{U(C)\\}{\rm tr}\\{U(C^{\prime})\\}$ $\displaystyle{\rm
for\leavevmode\nobreak\ }i=3,4,6,7$

where $\Gamma_{i}$ are all unique loops for a given shape, including loops in
the perpendicular direction – i.e. chair-shaped – for $i=$1,2,3,4.

Figure 2: Accuracy on the gradients in $\beta$ for the Wilson loops, $\mathcal{W}_{i}$ depicted in Fig. 1, measured using the perturbative expansion (PE) in Eq. (3) or alternatively using the new approach in Eq. (6) with the LO or the NLO terms for the flow action $\tilde{S}$ in Eq. (24). We depict the deviation $\Delta_{T}$ from the values measured using PE normalized such that the errors of PE are unitary, see Eq. (26). Thus, the continuous lines are the central values of PE, while the dashed lines are plus or minus one unit of the error of PE. Results at various $\beta$ are computed over ten thousand configurations with lattice size $16^{4}$, simulated with HMC and separated by 4 MDUs. |  |  $\beta=0$ |  $\beta=1$ |  $\beta=2$ |  $\beta=3$ |  $\beta=4$ |  $\beta=5$ |  $\beta=6$
---|---|---|---|---|---|---|---|---
LO | min | 0.08  % | 12.0  % | 23.1  % | 37.6  % | 45.4  % | 67.6  % | 65.3  %
max | 0.45  % | 16.7  % | 33.5  % | 48.2  % | 66.2  % | 81.1  % | 92.5  %
| $\\!\\!\\!\\!\\!\\!\sigma_{C}/\sigma_{S^{\prime}}$ | 0.00  % | 15.0  % | 29.2  % | 47.1  % | 65.6  % | 75.6  % | 102.0  %
NLO | min | 0.08  % | 1.8  % | 6.9  % | 17.0  % | 29.6  % | 55.3  % | 68.9  %
max | 0.45  % | 2.5  % | 10.4  % | 22.4  % | 43.9  % | 66.3  % | 99.1  %
| $\\!\\!\\!\\!\\!\\!\sigma_{C}/\sigma_{S^{\prime}}$ | 0.00  % | 2.4  % | 9.7  % | 22.8  % | 42.3  % | 62.2  % | 100.7  %
Table 1: Smallest and largest value for the ratio of the error obtained using
the new approach defined by Eq. (6) and the error of the perturbative
expansion given by Eq. (3) over the measured Wilson loops. Results are given
for the cases when the LO or the NLO is used in the flow action. We report
also the ratio between the standard deviation of $C$ and the one of
$S^{\prime}=\partial S/\partial\beta$, noting a strong connection between this
ratio and the one between the errors.

### Gradients of Wilson loops.

Since the convergence of the $\beta$-expansion in Eq. (24) deteriorates as
$\beta$ grows, we show results for simulations performed at various values of
$\beta$ from $\beta=0$ up to $\beta=6$. At $\beta=0$ the LO of the
$\beta$-expansion solves Eq. (23) exactly, while, on the other hand, the
region of physical interest starts from $\beta\gtrsim 5.8$, where a lattice
spacing of about $0.14$ fm is measured [19]. We compute the slope in $\beta$
for a set of Wilson loops, $\mathcal{W}_{i}$, using either the perturbative
expansion (PE) or the new approach using the LO or the NLO of Eq. (24)
computed at each $\beta$. For all experiments, we use a lattice of size
$16^{4}$ and measure expectation values over ten thousand configurations
simulated at each value of $\beta$ using an HMC algorithm. To compare
gradients for various Wilson loops, we study the deviation of the value
measured using an approach “$T$” from the weighted average of all approaches,
namely

$\Delta_{T,i}=\frac{1}{\sigma_{\rm PE,i}}\left(\frac{d\langle
W_{i}\rangle}{d\beta}\Bigg{|}_{T}-\frac{d\langle
W_{i}\rangle}{d\beta}\Bigg{|}_{\rm PE}\right)\,,$ (26)

where $\frac{d\langle W_{i}\rangle}{d\beta}\big{|}_{T}$ is the value from the
approach $T$ and we normalize the results using the error from the
perturbative expansion $\rm PE$, such that the error of $\Delta_{{\rm PE},i}$
is always unity. Results are depicted in Fig. 2 and summarized in Table 1. As
can be seen, we observe excellent statistical agreement between the results
measured using the perturbative expansion and the novel approach discussed in
this work, supporting the correctness of the derivation. We also observe an
impressive improvement in the error achieved at small $\beta$. This is
highlighted in Table 1, in which we list the minimum and maximum values of the
ratio between the errors obtained using the various approaches. At $\beta=0$,
results for LO and NLO are the same, since LO solves exactly Eq. (23). Here we
achieve the largest improvement for the signal-to-noise ratio, which is of
more than three orders of magnitude in calculating the slope in $\beta$ of the
plaquette. Additionally, looking at Table 1, we observe a very strong
connection between the improvement in the errors and the ratio of the standard
deviations of $C$ and $\partial S/\partial\beta$. The latter, therefore, can
be used to have an estimate of the gain one would achieve in a generic
situation.

Figure 3: Study of the convergence of the $\beta$-expansion, showing ratios
between the standard deviation of the LO, the NLO, and the next-to-next-
leading-order (NNLO) of the flow action or its Laplacian. Lines are
interpolations for guiding the eye.

### Convergence of the $\beta$-expansion.

On a side note, our experiments let us comment on the convergence of the
$\beta$-expansion in Eq. (24). This provides some interesting insights into
the ability to flow in $\beta$ using gradient flows [20]. On one side, if all
$\tilde{S}^{(k)}$ have similar magnitude, the series would be converging only
for $\beta<6$. On the other hand, $\tilde{S}^{(k)}$ might decrease in
magnitude at increasing $k$, therefore converging for all $\beta$. In Ref.
[1], an analytical study on the convergence of the $t$-expansion is found and
similar arguments can be applied to the $\beta$-expansion since they share the
same origin. In the following, on the other hand, we provide numerical
evidence for the first few orders. In our numerical experiments we have
computed the actions $\tilde{S}^{(0)}$ and $\tilde{S}^{(1)}$. Additionally, we
can measure the value of $\mathcal{L}_{0}\tilde{S}^{(2)}$ from its definition,
namely

$\mathcal{L}_{0}\tilde{S}^{(2)}=\left(\partial\mathcal{W}_{0},\partial\tilde{S}^{(1)}\right)\,,$
(27)

which is the scalar product of the forces of $\mathcal{W}_{0}$ and of
$\tilde{S}^{(1)}$, both computed in our experiments. Results for the ratio of
the first three orders are given in Fig. 3, in which we analyze the standard
deviation of these terms and not their value since the flow action is defined
up to a constant (23). While $\sigma(\tilde{S}^{(1)})$ is almost two orders of
magnitude smaller than $\sigma(\tilde{S}^{(0)})$, the ratio
$\sigma(\mathcal{L}_{0}\tilde{S}^{(2)})/\sigma(\mathcal{L}_{0}\tilde{S}^{(0)})$
is larger than
$\sigma(\mathcal{L}_{0}\tilde{S}^{(1)})/\sigma(\mathcal{L}_{0}\tilde{S}^{(0)})$
and increases in $\beta$, pointing in the direction that the various orders
might have similar magnitude. If this is true, then the $\beta$-expansion is
converging only for $\beta<6$ while at larger $\beta$ possibly all terms
contribute to the flow action, making it an extensive operator in the region
where the regime of physical interest begins. This is unfortunate for flow-
based trivializing maps since flowing at large $\beta$ becomes very
challenging. Indeed, to date, only flowing at very small $\beta$ has been
proven successful in four dimensions [21].

##  V Conclusions

In this work, we have presented a novel approach for computing gradients of
observables with respect to action parameters, which, in its simpler form
given in Eq. (8), is free of any disconnected contribution or vacuum
expectation value. We provide numerical evidence within our case study that
when a flow action solution to Eq. (9) is found, the new approach provides a
significant improvement in the signal-to-noise ratio compared to the commonly-
used perturbative expansion. We use, as a case study, the calculation of the
slope in $\beta$ of Wilson loops in pure-gauge SU(3) Yang-Mills theories. In
this situation, we are able to provide in Eq. (24) an analytic solution to the
required flow action in terms of a series expansion, referred to as
$\beta$-expansion. We compute the first two orders of the series and show
results at various $\beta$ for many Wilson loops. At $\beta=0$, the expansion
converges at the leading order, and in this case, we obtain the largest
improvement in the errors. At finite $\beta$, instead, the expansion converges
in orders of $\beta/6$. Therefore, a reduction in the errors is seen only for
$\beta<6$, which is a region of non-physical interest. On the other hand,
these results are only meant as proof of concept, since the approach will
change according to the application. Indeed, each case should be studied
independently, searching for an appropriate flow action, that satisfies Eq.
(9) for a desired action $S$ and parameter of the action $\theta$. Where
possible, this can be done analytically or, alternatively, it can be done
numerically using a parametric definition of $\tilde{S}$ [20] and machine
learning techniques to find the optimal flow action. In this case, for the
training, we suggest minimizing either the variance of
$\langle\mathcal{C}\rangle$, see Eq. (7), since the final objective is to make
it constant, i.e. Eq. (9); or to perform an observable-dependent training by
minimizing the value or the variance of
$\langle\mathcal{O}\rangle\langle\mathcal{C}\rangle-\langle\mathcal{O}\mathcal{C}\rangle$.
We have tested such an approach for the case study presented here and we have
obtained results compatible with the analytical solution and no further
improvement was observed. On the other hand, for all other cases where an
analytical solution is not easy to find, this is a new valuable application of
machine-learning techniques.

###### Acknowledgements.

### Acknowledgments.

The author wishes to thank Pan Kessel, Stefan Schaefer, and Lorenz Vaitl for
useful discussions and the enjoyable collaboration on gradient flows. The
author acknowledges financial support from the Cyprus Research and Innovation
Foundation under the project QC4LGT (grant agreement No. EXCELLENCE/0421/0019)
and 3D-NUCLEON (grant agreement No. EXCELLENCE/0421/0043). Numerical results
were obtained using the Cyclamen cluster of The Cyprus Institute equipped with
P100 GPUs. The software is implemented using the Lyncs-API [22], its interface
to QUDA [23] and the QUDA library [24]. All figures are available under a CC
BY 4.0 license111https://creativecommons.org/licenses/by/4.0/.

## References

  * Lüscher [2010] M. Lüscher, Trivializing maps, the Wilson flow and the HMC algorithm, Commun. Math. Phys. 293, 899 (2010), arXiv:0907.5491 [hep-lat] .
  * Horsley  _et al._ [2012] R. Horsley, R. Millo, Y. Nakamura, H. Perlt, D. Pleiter, P. Rakow, G. Schierholz, A. Schiller, F. Winter, and J. Zanotti, A lattice study of the glue in the nucleon, Physics Letters B 714, 312 (2012).
  * Chambers  _et al._ [2014] A. J. Chambers _et al._ (CSSM, QCDSF/UKQCD), Feynman-Hellmann approach to the spin structure of hadrons, Phys. Rev. D 90, 014510 (2014), arXiv:1405.3019 [hep-lat] .
  * Chambers  _et al._ [2015] A. J. Chambers _et al._ , Disconnected contributions to the spin of the nucleon, Phys. Rev. D 92, 114517 (2015), arXiv:1508.06856 [hep-lat] .
  * Bouchard  _et al._ [2017] C. Bouchard, C. C. Chang, T. Kurth, K. Orginos, and A. Walker-Loud, On the Feynman-Hellmann Theorem in Quantum Field Theory and the Calculation of Matrix Elements, Phys. Rev. D 96, 014504 (2017), arXiv:1612.06963 [hep-lat] .
  * Chang  _et al._ [2018] C. C. Chang _et al._ , A per-cent-level determination of the nucleon axial coupling from quantum chromodynamics, Nature 558, 91 (2018), arXiv:1805.12130 [hep-lat] .
  * Can  _et al._ [2020] K. U. Can _et al._ , Lattice QCD evaluation of the Compton amplitude employing the Feynman-Hellmann theorem, Phys. Rev. D 102, 114505 (2020), arXiv:2007.01523 [hep-lat] .
  * de Divitiis  _et al._ [2012] G. M. de Divitiis _et al._ , Isospin breaking effects due to the up-down mass difference in Lattice QCD, JHEP 04, 124, arXiv:1110.6294 [hep-lat] .
  * de Divitiis  _et al._ [2013] G. M. de Divitiis, R. Frezzotti, V. Lubicz, G. Martinelli, R. Petronzio, G. C. Rossi, F. Sanfilippo, S. Simula, and N. Tantalo (RM123), Leading isospin breaking effects on the lattice, Phys. Rev. D 87, 114505 (2013), arXiv:1303.4896 [hep-lat] .
  * Borsanyi  _et al._ [2015] S. Borsanyi _et al._ (BMW), Ab initio calculation of the neutron-proton mass difference, Science 347, 1452 (2015), arXiv:1406.4088 [hep-lat] .
  * Giusti  _et al._ [2017] D. Giusti, V. Lubicz, C. Tarantino, G. Martinelli, F. Sanfilippo, S. Simula, and N. Tantalo, Leading isospin-breaking corrections to pion, kaon and charmed-meson masses with Twisted-Mass fermions, Phys. Rev. D 95, 114504 (2017), arXiv:1704.06561 [hep-lat] .
  * Borsanyi  _et al._ [2021] S. Borsanyi _et al._ , Leading hadronic contribution to the muon magnetic moment from lattice QCD, Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * Giusti  _et al._ [2018] D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo, S. Simula, N. Tantalo, and C. Tarantino, First lattice calculation of the QED corrections to leptonic decay rates, Phys. Rev. Lett. 120, 072001 (2018), arXiv:1711.06537 [hep-lat] .
  * Di Carlo  _et al._ [2019] M. Di Carlo, D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo, S. Simula, and N. Tantalo, Light-meson leptonic decay rates in lattice QCD+QED, Phys. Rev. D 100, 034514 (2019), arXiv:1904.08731 [hep-lat] .
  * Alexandrou  _et al._ [2023] C. Alexandrou _et al._ (Extended Twisted Mass), Lattice calculation of the short and intermediate time-distance hadronic vacuum polarization contributions to the muon magnetic moment using twisted-mass fermions, Phys. Rev. D 107, 074506 (2023), arXiv:2206.15084 [hep-lat] .
  * Dragos  _et al._ [2021] J. Dragos, T. Luu, A. Shindler, J. de Vries, and A. Yousif, Confirming the Existence of the strong CP Problem in Lattice QCD with the Gradient Flow, Phys. Rev. C 103, 015202 (2021), arXiv:1902.03254 [hep-lat] .
  * Alexandrou  _et al._ [2021] C. Alexandrou, A. Athenodorou, K. Hadjiyiannakou, and A. Todaro, Neutron electric dipole moment using lattice QCD simulations at the physical point, Phys. Rev. D 103, 054501 (2021), arXiv:2011.01084 [hep-lat] .
  * Bhattacharya  _et al._ [2021] T. Bhattacharya, V. Cirigliano, R. Gupta, E. Mereghetti, and B. Yoon, Contribution of the QCD $\Theta$-term to the nucleon electric dipole moment, Phys. Rev. D 103, 114507 (2021), arXiv:2101.07230 [hep-lat] .
  * Schaefer  _et al._ [2011] S. Schaefer, R. Sommer, and F. Virotta (ALPHA), Critical slowing down and error analysis in lattice QCD simulations, Nucl. Phys. B 845, 93 (2011), arXiv:1009.5228 [hep-lat] .
  * Bacchio  _et al._ [2023] S. Bacchio, P. Kessel, S. Schaefer, and L. Vaitl, Learning trivializing gradient flows for lattice gauge theories, Phys. Rev. D 107, L051504 (2023), arXiv:2212.08469 [hep-lat] .
  * Abbott  _et al._ [2023] R. Abbott _et al._ , Normalizing flows for lattice gauge theory in arbitrary space-time dimension, (2023), arXiv:2305.02402 [hep-lat] .
  * Bacchio  _et al._ [2022] S. Bacchio, J. Finkenrath, and C. Stylianou, Lyncs-API: a Python API for Lattice QCD applications, PoS LATTICE2021, 542 (2022), arXiv:2201.03873 [hep-lat] .
  * Yamamoto  _et al._ [2022] S. Yamamoto, S. Bacchio, and J. Finenrath, Running HMC Simulation with Python via QUDA, (2022), arXiv:2212.06657 [hep-lat] .
  * Clark  _et al._ [2010] M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Solving Lattice QCD systems of equations using mixed precision solvers on GPUs, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .