#  Four Gluon Vertex from Lattice QCD

Manuel Colaço<EMAIL_ADDRESS>Orlando Oliveira<EMAIL_ADDRESS>Paulo J. Silva<EMAIL_ADDRESS>CFisUC, Departamento de Física, Faculdade de
Ciências e Tecnologia, Universidade de Coimbra, 3004-516 Coimbra, Portugal

###### Abstract

A lattice QCD calculation for the four gluon one-particle irreducible Green
function in the Landau gauge is discussed. Results for some of the associated
form factors are reported for kinematical configurations with a single
momentum scale. Our results show that the computation of this Green function
requires large statistical ensembles with 10K or larger number of gauge
configurations. The simulations considered herein have a clear Monte Carlo
signal for momenta up to $\sim 1$ GeV. The form factors show an hierarchy,
with the form factor associated with the tree level Feynman rule being
dominant and essentially constant for the range of momenta accessed. The
remaining form factors seem to increase as the momentum decreases, suggesting
that a possible $\log$ divergence may occur. The computed form factors are, at
least, in qualitative agreement with the results obtained with continuum
approaches to this vertex, when available.

###### Contents

  1. I Introduction
  2. II The Four Gluon Vertex
    1. II.1 An orthogonal set of operators
    2. II.2 Full Green function and lattice form factors
    3. II.3 Lattice form factors $F^{(i)}$ and 1-PI form factors
  3. III Lattice setup and choice of momenta
  4. IV Lattice form factors
  5. V Summary and Conclusions
  6. A Three gluon contributions for proportional momenta
  7. B Implementation of color-Lorentz projectors to measure $F^{(i)}$

##  I Introduction

The dynamical content of the QCD Green functions is of upmost importance to
all hadronic physics, specially its non-perturbative contents. Indeed, among
other phenomena, the interpretation of the hadron spectra, the comprehension
of dynamical chiral symmetry breaking and of the quark and gluon confinement
mechanisms are outside the scope of the perturbative solution of the theory.

The analysis of the QCD Green functions started looking at those functions
with smaller number of external legs, namely the propagator or two point
correlation functions. Their tensorial structure is simpler, as it requires
smaller tensor basis, and can be guided with the help of the Slavnov-Taylor
identities. All the two point QCD correlation functions have been studied in
detail with non-perturbative methods and a fairly good understanding of these
functions has been achieved Cucchieri and Mendes (2007); Aguilar  _et al._
(2008); Fischer  _et al._ (2009); Bogolubsky  _et al._ (2009); Oliveira and
Silva (2012); Ayala  _et al._ (2012); Duarte  _et al._ (2017). A good and
coherent picture of the two point QCD correlation functions, obtained with
different non-perturbative methods, has been achieved but there are issues
that are still under debate. An example is their analytical structure,
required in the computation of meson and baryon properties with continuum
methods, that remains to be understood Alkofer  _et al._ (2004); Strauss  _et
al._ (2012); Falcão  _et al._ (2020); Falcão and Oliveira (2022); Boito  _et
al._ (2023).

In QCD the three point one-particle irreducible (1-PI) Green functions have
also been object of study. Indeed, the ghost-gluon vertex Ilgenfritz  _et al._
(2007); Schleifenbaum  _et al._ (2005); Cucchieri  _et al._ (2008); Huber and
von Smekal (2013), the three gluon vertex Binosi  _et al._ (2013); Boucaud
_et al._ (2014); Eichmann  _et al._ (2014); Athenodorou  _et al._ (2016);
Duarte  _et al._ (2016a); Rodríguez-Quintero  _et al._ (2017); Binosi and
Papavassiliou (2018); Vujinovic and Mendes (2019); Aguilar  _et al._ (2019,
2020); Maas and Vujinović (2022); Aguilar  _et al._ (2021); Figueroa and
Peláez (2022); Catumba  _et al._ (2022); Barrios  _et al._ (2022); Aguilar
_et al._ (2023a); Pinto-Gómez  _et al._ (2023); Pinto-Gómez (2023); Aguilar
_et al._ (2023b) and the quark-gluon vertex Davydychev  _et al._ (2001);
Bender  _et al._ (2002); Skullerud and Kizilersu (2002); Skullerud  _et al._
(2003); Bhagwat and Tandy (2004); Kizilersu  _et al._ (2007); Matevosyan  _et
al._ (2007); Alkofer  _et al._ (2009); He (2009); Rojas  _et al._ (2013);
Aguilar  _et al._ (2014); Peláez  _et al._ (2015); Binosi  _et al._ (2017);
Bermudez  _et al._ (2017); Oliveira  _et al._ (2018, 2019); Sultan  _et al._
(2021); Kızılersü  _et al._ (2021); Gao  _et al._ (2021); Skullerud  _et al._
(2022); El-Bennich  _et al._ (2022); Marques  _et al._ (2023); Aguilar  _et
al._ (2023c) have been investigated with some degree of detail and its main
features are now understood. If for some of the vertices a qualitative picture
of their dependence for different kinematics can be claimed, a complete and
accurate description is still to be done.

The next level of QCD one-particle irreducible Green functions to be studied
with non-perturbative methods are those with four external legs. These include
the four quark and the four gluon correlation functions. The four quark Green
function is required to understand, for example, the meson spectra. The four
gluon correlation function has a similar role as the four quark function but
for the glueball spectra. Moreover, it also allows for the definition of a
renormalization group invariant QCD charge, it contributes to the Dyson-
Schwinger equation for the gluon propagator and impacts on the gluon
propagator itself Eichmann  _et al._ (2021). Herein, our focus of attention is
the four gluon one-particle irreducible Green function.

The non-perturbative content of the four gluon 1-PI Green function has been
studied with the continuum formulation of QCD using the Dyson-Schwinger
equations Kellermann and Fischer (2008); Binosi  _et al._ (2014); Cyrol  _et
al._ (2015). A full description of this four leg function requires a tensor
basis with a large number of elements Gracey (2014); Binosi  _et al._ (2014);
Eichmann  _et al._ (2015), that makes the analysis of this four point
correlation function difficult. Continuum studies of this function have to
make simplifications to arrive at a manageable computational scenario and the
calculations performed have to consider a basis with a reduced number of
elements. The continuum based Dyson-Schwinger calculations identify a number
of form factors that display very mild dependence with the momentum, as seems
to be the case of the form factor associated with the tree level tensor
structure that, nevertheless, seems to be suppressed at low momenta, and
others that show $\log$ divergences in the IR region. These divergences are
associated with the unprotected $\log$’s coming from the massless ghost loops
in the Dyson-Schwinger equations Binosi  _et al._ (2014). The continuum QCD
methods based on the Dyson-Schwinger equation for this 1-PI Green function
seem to favour such type of scenario.

The picture that emerges from the various continuum approaches to the four
gluon 1-PI Green function, although being consistent at the qualitative level,
also show differences that probably come from the different truncations and/or
the required modelling to solve the underlying equations. The differences
between the various studies provide an extra motivation to a lattice
calculation of the four gluon vertex. Indeed, the interplay of the various
non-perturbative methods is important to achieve a proper and sufficiently
accurate description of the QCD one-particle irreducible Green function.

In this work we address the computation of the four gluon one-particle
irreducible Green function within pure Yang-Mills lattice QCD simulations in
the Landau gauge. Preliminary results can be found in Colaço  _et al._ (2023).
The manuscript is organized as follows. In Sec. II it is given a description
of the four gluon full Green function, accessed in the lattice QCD
simulations, and how it relates to the corresponding 1-PI function.
Furthermore, details on the extraction of the 1-PI contribution in a
simulation in the Landau gauge are also discussed there, together with the
tensor basis that will be used for the calculation. The definitions of the
form factors to be measured can also be found in this section. The lattice
setup, the choice of the kinematic configurations to be used in the
measurements and some details of the lattice calculation are given in Sec.
III. The lattice form factors, including those associated with the amputated
Green functions, are reported and discussed in Sec. IV. Finally in Sec. V we
summarize and conclude. Some of the details of the calculation, namely a proof
of the decoupling of the three gluon contribution to the four point function
for a certain class of kinematical configuration and the lattice version of
the continuum operators, are given in the appendices.

##  II The Four Gluon Vertex

Let us start by discussing the computation of the four gluon 1-PI Green
function using the continuum formulation of QCD in Minkowski spacetime. By
doing so we aim to arrive at a better understanding of this complex Green
function before plugging into the measurement of the Green function with
lattice QCD simulations.

Figure 1: Diagramatic representation of the four gluon full Green function
defined in Eq (4). The filled vertex in the left is the complete four point
Green function, the blob filled in back is the connected Green function and
the blobs in grey stand for one particle-irreducible diagrams. Color, Lorentz
and momentum indices are omitted.

In the numerical Monte Carlo simulations of QCD on a lattice only the full
Green functions are accessed. In real space, these Green functions are defined
as the vacuum expectation values of the time order product of the fundamental
fields. In QCD, the gluon N-point full Green function reads

$\mathcal{G}^{a_{1}\cdots
a_{N}}_{\mu_{1}\cdots\mu_{N}}(x_{1},\,\cdots,\,x_{N})=\langle~{}A^{a_{1}}_{\mu_{1}}(x_{1})~{}A^{a_{2}}_{\mu_{2}}(x_{2})~{}\cdots~{}A^{a_{N}}_{\mu_{N}}(x_{N})~{}\rangle\
,$ (1)

where $\langle\cdots\rangle$ represent the vaccuum expectation value, realised
as an ensemble average over gauge configurations sampled with an appropriate
action. Any other quantity, like form factors, have to be computed a
posteriori by handling a set of full Green functions $\mathcal{G}^{a_{1}\cdots
a_{N}}_{\mu_{1}\cdots\mu_{N}}$. The computation of the 1-PI four gluon Green
function is no exception and the first step to perform is to decompose the
4-point full Green function in terms of QCD 1-PI Green functions. For pure
Yang-Mills theory the functional generator of the full Green functions $Z[J]$
is

$Z[J]=e^{i\,W[J]}=\int\mathcal{D}A\,e^{i\,S[A]+i(J,A)}\ ,$ (2)

where $S[A]$ is the effective action associated with the Yang-Mills theory,
that has contributions coming from the gauge and ghost fields,

$(J,A)=\int d^{4}x~{}A^{a}_{\mu}(x)J^{a\,\mu}(x)\ ,$ (3)

$W[J]$ is the generating functional of the connected Green functions and $J$
are classical external sources. It follows from the definition of the full
Green function that the four gluon full Green function is given by

$\displaystyle\mathcal{G}^{a_{1}\cdots
a_{4}}_{\mu_{1}\cdots\mu_{4}}(x_{1},\,\cdots,\,x_{4})$ $\displaystyle=$
$\displaystyle\left.\frac{\delta^{4}Z[J]}{i\,\delta
J^{a_{1}\,\mu_{1}}(x_{1})~{}\cdots~{}i\,\delta
J^{a_{4}\,\mu_{4}}(x_{4})}\right|_{J=0}$ (4) $\displaystyle=$
$\displaystyle\left.i\,\frac{\delta^{4}W[J]}{\delta
J^{a_{1}\,\mu_{1}}(x_{1})~{}\cdots~{}\delta
J^{a_{4}\,\mu_{4}}(x_{4})}\right|_{J=0}$
$\displaystyle~{}-\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{1}\,\mu_{1}}(x_{1})~{}\delta
J^{a_{2}\,\mu_{2}}(x_{2})}\right|_{J=0}~{}\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{3}\,\mu_{3}}(x_{3})~{}\delta J^{a_{4}\,\mu_{4}}(x_{4})}\right|_{J=0}$
$\displaystyle\quad\quad-\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{1}\,\mu_{1}}(x_{1})~{}\delta
J^{a_{3}\,\mu_{3}}(x_{3})}\right|_{J=0}~{}\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{2}\,\mu_{2}}(x_{2})~{}\delta J^{a_{4}\,\mu_{4}}(x_{4})}\right|_{J=0}$
$\displaystyle\qquad\qquad-\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{1}\,\mu_{1}}(x_{1})~{}\delta
J^{a_{4}\,\mu_{4}}(x_{4})}\right|_{J=0}~{}\left.\frac{\delta^{2}W[J]}{\delta
J^{a_{2}\,\mu_{2}}(x_{2})~{}\delta J^{a_{3}\,\mu_{3}}(x_{3})}\right|_{J=0}\ .$

The diagrammatic representation of $\mathcal{G}^{a_{1}\cdots
a_{4}}_{\mu_{1}\cdots\mu_{4}}(x_{1},\,\cdots,\,x_{4})$ in terms of connected
Green functions is given in Fig. 1. The four gluon full Green function is a
sum of a four-gluon connected Green function, represented by the full black
blob in Fig. 1, and disconnected parts that are functions of the two gluon
Green function, i.e. of the gluon propagator, the full blobs in grey. The four
gluon connected Green function itself can be written in terms of 1-PI
functions. In order to workout this decomposition let us consider the
classical Yang-Mills field given by

$A_{cl,\mu}^{a}(x)=\frac{\delta W[J]}{\delta J^{a\,\mu}(x)}\ $ (5)

and introduce the generating functional for the one-particle irreducible
correlation functions

$\Gamma[A_{cl}]=W[J]-(J,A_{cl})\ .$ (6)

The one-particle irreducible Green functions are the functional derivatives of
the $\Gamma[A_{cl}]$ at vanishing classical fields. It follows from the above
definitions that

$\frac{\delta\,\Gamma}{\delta A_{cl,\mu}^{a}(x)}=-\,J^{a\,\mu}(x)$ (7)

and, from this relation combined with Eqs (5) and (6), the orthogonality
relation

$\int d^{4}z~{}~{}\frac{\delta^{2}W[J]}{\delta J^{a}_{\mu}(x)~{}~{}\delta
J^{c}_{\nu}(z)}~{}~{}\frac{\delta^{2}\,\Gamma[A_{cl}]}{\delta
A_{cl}^{c}\,{}^{\nu}(z)~{}~{}\delta
A_{cl}^{b}\,{}_{\zeta}(y)}=-\,\delta^{ab}~{}g^{\mu\zeta}~{}\delta(x-y)$ (8)

is derived. By taking functional derivatives of Eq (7), using the
orthogonality relation (8), and after some straightforward algebra, one
arrives at the equality

$\displaystyle\left.\frac{\delta^{4}W[J]}{\delta J^{a}_{\mu}(x)~{}\delta
J^{b}_{\nu}(y)~{}\delta J^{c}_{\iota}(z)~{}\delta
J^{e}_{\zeta}(t)}\right|_{J=0}~{}=~{}\int
d^{4}w_{1}~{}d^{4}w_{2}~{}d^{4}w_{3}~{}\Bigg{\\{}$ $\displaystyle\qquad
D^{bb^{\prime}}_{\nu\nu^{\prime}}(y-w_{1})~{}D^{ee^{\prime}}_{\zeta\zeta^{\prime}}(t-w_{2})~{}\frac{\delta^{3}W}{\delta
J^{a}_{\mu}(x)~{}\delta J^{c}_{\iota}(z)~{}\delta
J^{c^{\prime}}_{\iota^{\prime}}(w_{3})}~{}\frac{\delta^{3}\Gamma}{\delta
A^{b^{\prime}}_{cl\,\nu^{\prime}}(w_{1})~{}\delta
A^{c^{\prime}}_{cl\,\iota^{\prime}}(w_{3})~{}\delta
A^{e^{\prime}}_{cl\,\zeta^{\prime}}(w_{2})}$
$\displaystyle\qquad+~{}D^{cc^{\prime}}_{\iota\iota^{\prime}}(z-w_{1})~{}D^{ee^{\prime}}_{\zeta\zeta^{\prime}}(t-w_{2})~{}\frac{\delta^{3}W}{\delta
J^{a}_{\mu}(x)~{}\delta J^{b}_{\nu}(y)~{}\delta
J^{b^{\prime}}_{\nu^{\prime}}(w_{3})}~{}\frac{\delta^{3}\Gamma}{\delta
A^{b^{\prime}}_{cl\,\nu^{\prime}}(w_{3})~{}\delta
A^{c^{\prime}}_{cl\,\iota^{\prime}}(w_{1})~{}\delta
A^{e^{\prime}}_{cl\,\zeta^{\prime}}(w_{2})}$
$\displaystyle\qquad+~{}D^{cc^{\prime}}_{\iota\iota^{\prime}}(z-w_{1})~{}D^{ee^{\prime}}_{\nu\nu^{\prime}}(y-w_{2})~{}\frac{\delta^{3}W}{\delta
J^{a}_{\mu}(x)~{}\delta J^{b}_{\zeta}(t)~{}\delta
J^{b^{\prime}}_{\zeta^{\prime}}(w_{3})}~{}\frac{\delta^{3}\Gamma}{\delta
A^{b^{\prime}}_{cl\,\zeta^{\prime}}(w_{3})~{}\delta
A^{c^{\prime}}_{cl\,\iota^{\prime}}(w_{1})~{}\delta
A^{e^{\prime}}_{cl\,\nu^{\prime}}(w_{2})}~{}\Bigg{\\}}$ $\displaystyle+\int
d^{4}w_{1}~{}d^{4}w_{2}~{}d^{4}w_{3}~{}d^{4}w_{4}~{}D^{aa^{\prime}}_{\mu\mu^{\prime}}(x-w_{1})~{}D^{bb^{\prime}}_{\nu\nu^{\prime}}(y-w_{2})~{}D^{cc^{\prime}}_{\iota\iota^{\prime}}(z-w_{3})~{}D^{ee^{\prime}}_{\zeta\zeta^{\prime}}(t-w_{4})\times$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\frac{\delta^{4}\Gamma}{\delta
A^{a^{\prime}}_{cl\,\mu^{\prime}}(w_{1})~{}\delta
A^{b^{\prime}}_{cl\,\nu^{\prime}}(w_{2})~{}\delta
A^{c^{\prime}}_{cl\,\iota^{\prime}}(w_{3})~{}\delta
A^{e^{\prime}}_{cl\,\zeta^{\prime}}(w_{4})}$ (9)

whose diagrammatic representation is given in Fig. 2. The connected four gluon
Green function is a sum of terms that include the four gluon 1-PI Green
function and contributions proportional to the three gluon 1-PI Green
functions.

Figure 2: Diagramatic representation of the four gluon connected Green
function (full black blob) appearing in Fig. 1 in terms of one particle-
irreducible functions (full grey blobs).

On a lattice simulation it is the full Green function
$\mathcal{G}^{a_{1}\cdots a_{4}}_{\mu_{1}\cdots\mu_{4}}$, see Eq. (4) and Fig.
1, that is measured. Its decomposition in terms of connected Green functions
calls for the four-gluon connected Green function and disconnected terms that
are proportional to products of gluon propagators. In the lattice formulation
of QCD, the Green functions are proportional to the lattice volume and,
therefore, the terms represented in Fig. 1 are such that the four gluon
connected part is proportional to the lattice volume $V$, while the terms
proportional to the propagators are proportional to $V^{2}$. The lattice
volume dependence of the terms that define $\mathcal{G}^{a_{1}\cdots
a_{4}}_{\mu_{1}\cdots\mu_{4}}$ make it almost impossible to have a good Monte
Carlo signal for the connected four gluon diagram, unless the contribution of
the disconnected parts vanish. This can be achieved considering only those
kinematical configurations such that the momenta $p_{i}$ and $p_{j}$
associated, respectively, with the external legs $i$ and $j$ are such that
$p_{i}+p_{j}\neq 0$ for all $i$ and $j$. From now on, it will be assumed that
these conditions on the momenta of the external legs are fulfilled and,
therefore, the contributions due to the disconnected diagrams can be
discarded.

The four gluon connected Green function, represented diagrammatically in Fig.
2, is given by a term that includes the four gluon 1-PI Green function and
diagrams that include three gluon 1-PI Green functions and gluon propagators.
Then, to access the four gluon 1-PI Green function from the connected four
gluon Green function, the three gluon contributions needs to be known.
However, such information on the three gluon contributions is currently not
available. One can either introduce some modelling for the three gluon 1-PI
Green function, in such a way that its contribution can be removed, or further
constrain the kinematical configuration of the external legs to suppress the
contributions that call for three-gluon one-particle irreducible diagrams. For
example, as discussed in Binosi  _et al._ (2014), this is case when the
incoming momenta are $p_{1}=p$, $p_{2}=p$, $p_{3}=p$ and $p_{4}=-3\,p$. A
first and preliminary analysis of this kinematical configuration using lattice
QCD simulations can be found in Catumba (2021).

In general, for the Landau gauge, due to the orthogonality of the gluon
propagator, for the kinematical configurations where the momenta of the
external legs are all proportional, i.e. for the class of momenta such that

$p_{1}=p,\qquad p_{2}=\eta\,p,\qquad p_{3}=\lambda\,p,\qquad
p_{4}=-(1+\eta+\lambda)\,p\ ,$ (10)

where $\eta$ and $\lambda$ are real numbers, the terms in the decomposition of
the four gluon connected Green function that require the contribution of the
three-gluon vertex vanish. The proof of this statement can be found in App. A.
Herein, in order to avoid any type of extra modelling, we choose to consider a
subset of the kinematical configurations that belong to the class of momenta
defined in Eq. (10). Although the modelling of the three gluon 1-PI is
avoided, our choice based on the kinematics constraint the type of information
that can be accessed in the simulation.

The choice of momenta is important to access the four gluon 1-PI Green
function without extra assumptions. However, the color-Lorentz structure of
this Green function is complex. For a general kinematical configuration, the
number of tensors required to define a tensor basis to fully describe the four
gluon 1-PI is large. For example, for the symmetric point, the four gluon
vertex requires more than one hundred different tensors for its full
description Gracey (2014). General discussions of the tensor basis for this
Green function can be found in Binosi  _et al._ (2014); Eichmann  _et al._
(2015). Ideally, one would access each of the form factors that multiply each
tensor of the complete basis for the Green function with the help of a
suitable projection operator. However, not only the lattice simulations
convolute the 1-PI functions with the gluon propagators but also the quality
of the signal-to-noise ratio in a Monte Carlo simulation favor certain types
of projecting operators. In principle, the signal-to-noise ratio can always be
improved by generating larger ensembles of gauge configurations, but a
practical computation has always a limited statistical ensemble.

A complete tensor basis describing the four gluon 1-PI includes several types
of operators that are proportional to the metric tensor and/or momenta.
However, given that in a lattice simulation the measured function convolutes
the four gluon 1-PI with gluon propagators, for the class of kinematical
configurations considered here that is characterized by a single momentum, see
Eq. (10), due to the orthogonality of the gluon propagator in the Landau
gauge, the tensors proportional to momentum do not contribute to the full
Green function that is measured in a simulation. This property simplifies
considerably the tensor analysis the Green function and the tensors that
contribute can only include the metric tensor for its Lorentz component and
$f_{abc}$ or $d_{abc}$, the fully anti-symmetric and fully symmetric structure
constants, and/or $\delta^{ab}$ in the description of the color part of the
Green function. Then, by imposing Bose symmetry it is possible to build the
possible tensor structures that, for the chosen kinematics, contribute to the
Green measured on a simulation. A particular tensor that belongs to the class
of allowed operators is the tensor structure that appears in the perturbative
tree level four gluon Feynman rule

$\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}=f_{abr}f_{cdr}\left(g_{\mu\eta}g_{\nu\zeta}-g_{\mu\zeta}g_{\nu\eta}\right)+f_{acr}f_{bdr}\left(g_{\mu\nu}g_{\eta\zeta}-g_{\mu\zeta}g_{\nu\eta}\right)+f_{adr}f_{bcr}\left(g_{\mu\nu}g_{\eta\zeta}-g_{\mu\eta}g_{\nu\zeta}\right)\
.$ (11)

From this tensor another operator that also contributes to the lattice Green
function can be build replacing the fully anti-symmetric SU(3) structure
constants $f_{abc}$ by the fully symmetric $d_{abc}$ and symmetrizing the
Lorentz structure part. The corresponding tensor reads

$\widetilde{\Gamma}^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}=d_{abr}d_{cdr}\left(g_{\mu\eta}g_{\nu\zeta}+g_{\mu\zeta}g_{\nu\eta}\right)+d_{acr}d_{bdr}\left(g_{\mu\zeta}g_{\nu\eta}+g_{\mu\nu}g_{\eta\zeta}\right)+d_{adr}d_{bcr}\left(g_{\mu\nu}g_{\eta\zeta}+g_{\mu\eta}g_{\nu\zeta}\right)\
,$ (12)

Besides the two above tensors, one may also consider the tensor operator

$\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}=\Big{(}\delta^{ab}\,\delta^{cd}+\delta^{ac}\,\delta^{bd}+\delta^{ad}\,\delta^{bc}\Big{)}~{}\Big{(}g_{\mu\nu}\,g_{\eta\zeta}+g_{\mu\eta}\,g_{\nu\zeta}+g_{\mu\zeta}\,g_{\nu\eta}\Big{)}$
(13)

that was mentioned in Binosi  _et al._ (2014). The three above tensors are not
orthogonal to each other and their color-Lorentz scalars are

$\displaystyle\widetilde{\Gamma}^{(0)}\,\cdot\,\widetilde{\Gamma}^{(0)}$
$\displaystyle=$ $\displaystyle 108\,N^{2}\,\Big{(}N^{2}-1\Big{)}=7776\ ,$
(14) $\displaystyle\widetilde{\Gamma}^{(0)}\,\cdot\,\widetilde{\Gamma}^{(1)}$
$\displaystyle=$ $\displaystyle-36\,\Big{(}N^{4}-5\,N^{2}+4\Big{)}=-1440\ ,$
(15) $\displaystyle\widetilde{\Gamma}^{(0)}\,\cdot\,\widetilde{\Gamma}^{(2)}$
$\displaystyle=$ $\displaystyle 0$ (16)
$\displaystyle\widetilde{\Gamma}^{(1)}\,\cdot\,\widetilde{\Gamma}^{(1)}$
$\displaystyle=$ $\displaystyle
12\left(17\,N^{4}\,-\,209\,N^{2}\,-\,\frac{496}{N^{2}}\,+\,688\right)=\frac{4640}{3}\
,$ (17)
$\displaystyle\widetilde{\Gamma}^{(1)}\,\cdot\,\widetilde{\Gamma}^{(2)}$
$\displaystyle=$
$\displaystyle\frac{288\left(N^{4}\,-\,5\,N^{2}\,+\,4\right)}{N}=3840\ ,$ (18)
$\displaystyle\widetilde{\Gamma}^{(2)}\,\cdot\,\widetilde{\Gamma}^{(2)}$
$\displaystyle=$ $\displaystyle 216\Big{(}N^{4}\,-\,1\Big{)}=17280\ ,$ (19)

where the numbers in the r.h.s. refer to the particular case where $N=3$. In
general, the above set of tensors do not define a complete tensor basis for
the Green function that is measured on a simulation and others should be
considered. However, in the current work only their contribution will be
measured.

###  II.1 An orthogonal set of operators

Although the operators $\widetilde{\Gamma}^{(0)}$, $\widetilde{\Gamma}^{(1)}$
and $\widetilde{\Gamma}^{(2)}$ are not orthogonal in the color-Lorentz space,
by a linear combination of these operators it is possible to build three
orthogonal operators from $\widetilde{\Gamma}^{(i)}$. Indeed, a
straightforward calculation shows that the tensor operators

$\displaystyle t^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}$ $\displaystyle=$
$\displaystyle\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}\ ,$ (20)
$\displaystyle t^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}$ $\displaystyle=$
$\displaystyle\frac{N^{2}-4}{3\,N^{2}}~{}~{}\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}+~{}\widetilde{\Gamma}^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}-~{}\frac{4\left(N^{2}-4\right)}{3\,N\,\left(N^{2}+1\right)}~{}~{}\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}$
(21) $\displaystyle=$
$\displaystyle\frac{5}{27}~{}~{}\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}+~{}\widetilde{\Gamma}^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}-~{}\frac{2}{9}~{}~{}\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}\
,$ $\displaystyle t^{(3)}\,{}^{abcd}_{\mu\nu\eta\zeta}$ $\displaystyle=$
$\displaystyle\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}$ (22)

are orthogonal in the color-Lorentz space and, therefore,

$t^{(i)}\cdot t^{(j)}~{}=~{}\mathcal{N}_{j}~{}~{}\delta^{ij}\ ,$ (23)

where the normalization factors $\mathcal{N}_{j}$ read

$\displaystyle\mathcal{N}_{1}$ $\displaystyle=$ $\displaystyle
108\,N^{2}\,\Big{(}N^{2}-1\Big{)}=7776\ ,$ (24) $\displaystyle\mathcal{N}_{2}$
$\displaystyle=$ $\displaystyle\frac{1280}{3}\ ,$ (25)
$\displaystyle\mathcal{N}_{3}$ $\displaystyle=$ $\displaystyle
216\Big{(}N^{4}\,-\,1\Big{)}=17280\ ;$ (26)

once more the numbers in the r.h.s. are for $N=3$. For a general kinematics,
the three tensors $t^{(1)}$, $t^{(2)}$ and $t^{(3)}$ do not define a complete
basis neither for the full Green function nor for the one-particle function.
The exception being the configuration of momenta such that
$p_{1}=p_{2}=p_{3}=p$ and $p_{4}=-3\,p$ where $t^{(1)}$, $t^{(2)}$ and
$t^{(3)}$ describe completely the correlation function Binosi  _et al._
(2014).

###  II.2 Full Green function and lattice form factors

For the special kinematics under discussion herein, the full Green functions
measured in lattice simulations gets only contributions from the connected
four gluon full Green function. In the Landau gauge, the full Green function
can be written as

$\displaystyle\mathcal{G}^{abcd}_{\mu\nu\eta\zeta}(p_{1},\,p_{2},\,p_{3},\,p_{4})$
$\displaystyle=$
$\displaystyle\Big{(}P_{\perp}(p_{1})\Big{)}_{\mu\mu^{\prime}}\,\Big{(}P_{\perp}(p_{2})\Big{)}_{\nu\nu^{\prime}}\,\Big{(}P_{\perp}(p_{3})\Big{)}_{\eta\eta^{\prime}}\,\Big{(}P_{\perp}(p_{4})\Big{)}_{\zeta\zeta^{\prime}}~{}D(p^{2}_{1})\,D(p^{2}_{2})\,D(p^{2}_{3})\,D(p^{2}_{4})$
(28)
$\displaystyle~{}\Bigg{(}F(p^{2}_{1},\dots)\,t^{(1)}\,{}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}~{}+~{}G(p^{2}_{1},\dots)\,t^{(2)}\,{}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}~{}+~{}H(p^{2}_{1},\dots)\,t^{(3)}\,{}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}~{}+~{}\cdots\Bigg{)}$

where $F(p^{2}_{1},\dots)$, $G(p^{2}_{1},\dots)$, $H(p^{2}_{1},\dots)$, etc,
are Lorentz scalar form factors, $\cdots$ represent the contribution of the
remaining components of the tensor basis, supposedly but not necessarilly
orthogonal to the space spanned by $t^{(1)}$ to $t^{(3)}$, and

$\Big{(}P_{\perp}(p)\Big{)}_{\mu\nu}=g_{\mu\nu}-\frac{p_{\mu}\,p_{\nu}}{p^{2}}$
(29)

is the orthogonal projector in momentum space that appears in the definition
of the Landau gauge gluon propagator

$D^{ab}_{\mu\nu}(p)=\delta^{ab}~{}\Big{(}P_{\perp}(p)\Big{)}_{\mu\nu}~{}D(p^{2})\
.$ (30)

The 1-PI form factors to be measured from the full Green function are

$\displaystyle F^{(0)}$ $\displaystyle=$
$\displaystyle-\,\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}~{}\mathcal{G}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}(p_{1},\,p_{2},\,p_{3},\,p_{4})~{}~{}g^{\mu\mu^{\prime}}\,g^{\nu\nu^{\prime}}\,g^{\eta\eta^{\prime}}\,g^{\zeta\zeta^{\prime}}\
,$ (31) $\displaystyle F^{(1)}$ $\displaystyle=$
$\displaystyle-\,\widetilde{\Gamma}^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}~{}\mathcal{G}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}(p_{1},\,p_{2},\,p_{3},\,p_{4})~{}~{}g^{\mu\mu^{\prime}}\,g^{\nu\nu^{\prime}}\,g^{\eta\eta^{\prime}}\,g^{\zeta\zeta^{\prime}}\
,\ $ (32) $\displaystyle F^{(2)}$ $\displaystyle=$
$\displaystyle-\,\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}~{}~{}\mathcal{G}^{abcd}_{\mu^{\prime}\nu^{\prime}\eta^{\prime}\zeta^{\prime}}(p_{1},\,p_{2},\,p_{3},\,p_{4})~{}~{}g^{\mu\mu^{\prime}}\,g^{\nu\nu^{\prime}}\,g^{\eta\eta^{\prime}}\,g^{\zeta\zeta^{\prime}}\
.$ (33)

In general, the $F^{(i)}$ are linear combinations of the form factors $F$,
$G$, $H$ $\cdots$ that describe the 1-PI four gluon function. It follows from
the definitions (31) – (33) that their r.h.s. are symmetric under interchange
of any pair of momenta and Bose symmetry demands that the $F^{(i)}$ can only
depend on the momenta through the combinations

$F^{(i)}\equiv
F^{(i)}\Big{(}p^{2}_{1}+p^{2}_{2}+p^{2}_{3}+p^{2}_{4}~{},~{}(p_{1}p_{2})+(p_{1}p_{3})+(p_{1}p_{4})+(p_{2}p_{3})+(p_{2}p_{4})+(p_{3}p_{4})\Big{)}\
.$ (34)

For the kinematical configurations investigated in the current work, where all
momenta are proportional to each other and there is a unique momentum scale,
one can write that $F^{(i)}\equiv F^{(i)}(p^{2})$ to simplify the notation.

###  II.3 Lattice form factors $F^{(i)}$ and 1-PI form factors

Our primer concern in this work is to compute the $F^{(i)}$. These functions
are not the form factors that describe the 1-PI four gluon vertex, see the
definition in Eq. (28), but are given by linear combinations of $F$, $G$, $H$,
$\dots$ Assuming that the functions $F$, $G$ and $H$ give the dominant
contributions to the form factors lattice form factors $F^{(i)}$, then for the
kinematical configurations
$(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(0,~{}~{}p,~{}~{}p,~{}~{}-2p)$
and
$(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(0,~{}~{}p,~{}~{}2p,~{}~{}-3p)$,
a straightforward calculation gives

$\displaystyle F^{(0)}(p^{2})$ $\displaystyle=$ $\displaystyle 7776\,F(p^{2})\
,$ (35) $\displaystyle F^{(1)}(p^{2})$ $\displaystyle=$
$\displaystyle-\,1440\,F(p^{2})~{}+~{}\frac{1280}{3}\,G(p^{2})~{}+~{}3840\,H(p^{2})\
,$ (36) $\displaystyle F^{(2)}(p^{2})$ $\displaystyle=$ $\displaystyle
17280\,H(p^{2})$ (37)

while for
$(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(p,~{}~{}p,~{}~{}p,~{}~{}-3p)$
one has

$\displaystyle F^{(0)}(p^{2})$ $\displaystyle=$ $\displaystyle 3888\,F(p^{2})\
,$ (38) $\displaystyle F^{(1)}(p^{2})$ $\displaystyle=$
$\displaystyle-\,720\,F(p^{2})~{}+~{}\frac{640}{3}\,G(p^{2})~{}+~{}2400\,H(p^{2})\
,$ (39) $\displaystyle F^{(2)}(p^{2})$ $\displaystyle=$ $\displaystyle
10800\,H(p^{2})\ .$ (40)

If $F$ and $H$ are directly related to $F^{(0)}$ and $F^{(2)}$, the form
factor $G(p^{2})$ can be accessed by the linear combination of all the lattice
form factors

$\frac{5}{27}\,F^{(0)}(p^{2})+F^{(1)}(p^{2})-\frac{2}{9}\,F^{(2)}(p^{2})=\left\\{\begin{array}[]{l@{\hspace{0.6cm}}l}\frac{1280}{3}\,G(p^{2})\
,\hfil\hskip 17.07182pt&\mbox{ all kinematics except
}(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(p,\,p,\,p,\,-3p)\\\ \hfil\hskip
17.07182pt&\\\ \frac{640}{3}\,G(p^{2})\ ,\hfil\hskip 17.07182pt&\mbox{ for
}(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(p,\,p,\,p,\,-3p)\
.\end{array}\right.$ (41)

We close this section by recalling the reader that for the kinematics
$(p_{1},~{}~{}p_{2},~{}~{}p_{3},~{}~{}p_{4})=(p,\,p,\,p,\,-3p)$ the three
tensor structures considered define a complete basis Binosi  _et al._ (2014).

##  III Lattice setup and choice of momenta

The main goal of this work is to measure the form factors that describe the
four gluon one-particle irreducible correlation function using lattice
simulations. The results described below refer to simulations that use the
Wilson action for $\beta=6.0$. For this bare coupling constant the
corresponding lattice spacing, measured from the string tension, is
$a=0.1016(25)$ fm or, equivalently, $1/a=1.943(48)$ GeV. In the conversion
from lattice to physical units we use the central values just quoted111Note
that there is an uncertainty of $\sim$2.5% in the scale setting towards
physical momentum. However, in the following, the uncertainty associated with
$1/a$ will be ignored and only the central value will be used. See Boucaud
_et al._ (2017); Duarte  _et al._ (2017) for discussions.. The computer
simulations were performed using the Chroma Edwards and Joo (2005) and PFFT
Pippig (2013) libraries and run on the Navigator supercomputer at the
University of Coimbra.

The form factors $F^{(i)}$ defined in Eqs (31) to (33) were computed with
several ensembles of gauge configurations. Our first try was to consider the
ensembles used in Oliveira and Silva (2012); Duarte  _et al._ (2016b); Dudal
_et al._ (2018) that have 2000 gauge configurations for the $64^{4}$ lattice
and 1801 gauge configurations for the $80^{4}$ lattice, both rotated to the
Landau gauge222 The details on the sampling and the rotation of the gauge
configurations towards the Landau gauge can be found in the above cited
works.. For these ensembles of configurations the signal-to-noise ratios are
quite poor, preventing from achieving any reliable measurement for the
$F^{(i)}$. To overcome the problem of the signal-to-noise ratio we have
generated 9038 gauge configurations, in the Landau gauge, for a $32^{4}$
lattice with the same $\beta$ value, together with 4560 gauge configurations
for a $48^{4}$ lattice and $\beta=6.0$. Our choice for $\beta$ relies on the
study of the gluon propagator Oliveira and Silva (2012) that suggests that for
this $\beta$ the lattice spacing effects are under control with possible mild
volume effects in the IR that show up only in the smallest lattice. From now
on, we will refer only to the measurements with these two latter lattices.

Our choice of momenta to measure the form factors $F^{(i)}$ considers only a
subset of the class of momenta defined in (10). The subset of momenta was
chosen with the aim to minimize the effects due to the breaking of rotational
symmetry and our choice was to use momenta $p_{i}$ as close to each other as
possible. The momentum configurations considered are

$\begin{array}[]{l@{\hspace{1cm}}l@{\hspace{1cm}}l@{\hspace{1cm}}l@{\hspace{1cm}}l@{\hspace{1cm}}}p_{1}=0,\hfil\hskip
28.45274pt&p_{2}=p,\hfil\hskip 28.45274pt&p_{3}=p,\hfil\hskip
28.45274pt&p_{4}=-2\,p\hfil\hskip 28.45274pt&~{}\mbox{ referred as
}~{}(0,~{}~{}p,~{}~{}~{}\,p,~{}~{}-2p)\ ,\hfil\hskip 28.45274pt\\\
p_{1}=0,\hfil\hskip 28.45274pt&p_{2}=p,\hfil\hskip
28.45274pt&p_{3}=2p,\hfil\hskip 28.45274pt&p_{4}=-3\,p\hfil\hskip
28.45274pt&~{}\mbox{ referred as }~{}(0,~{}~{}p,~{}~{}2p,~{}~{}-3p)\
,\hfil\hskip 28.45274pt\\\ p_{1}=p,\hfil\hskip 28.45274pt&p_{2}=p,\hfil\hskip
28.45274pt&p_{3}=p,\hfil\hskip 28.45274pt&p_{4}=-3\,p\hfil\hskip
28.45274pt&~{}\mbox{ referred as }~{}(p,~{}~{}p,~{}~{}~{}\,p,~{}~{}-3p)\
.\hfil\hskip 28.45274pt\end{array}$ (42)

In a simulation that uses hypercubic lattices, the breaking of rotational
symmetry implies that the measured form factors $F^{(i)}$ are not only
functions of $p^{2}$ but also of the H4 invariants $p^{[n]}$ Becirevic  _et
al._ (2000); de Soto and Roiesnel (2007); Vujinovic and Mendes (2019); Catumba
_et al._ (2021). If for the gluon propagator the building of these invariants
and the evaluation of the effects due to the breaking of rotational symmetry
is relatively straightforward, a similar analysis of Green functions with
higher number of external legs becomes rather involved. For the three point
functions, the effects due to the breaking of rotational symmetry were
discussed in Vujinovic and Mendes (2019) and for an hypercubic lattice, the
breaking of rotational symmetry is minimized when one of the momenta is
proportional to $(1,\,1,\,1,\,1)$. Similar conclusions have been observed for
the two point functions Leinweber  _et al._ (1999); de Soto and Roiesnel
(2007); Catumba  _et al._ (2021) At least for the two point and three point
gluon functions it was observed that the replacement of the naive momentum by
an improved momentum, see Eq. (43) in Sec. IV, seem to handle the effects due
to the breaking of rotation symmetry.

For the kinematical configurations under consideration, the form factors
$F^{(i)}$ are functions of a single momentum scale. In order to improve the
signal-to-noise ratio, given that the form factors are functions of $p^{2}$,
in all cases we average over equivalent momenta, including negative momenta,
before performing the ensemble averages. For example, for each gauge
configuration the form factors associated with the momenta

$\begin{array}[]{lrrrrr@{\hspace{0.4cm}}lrrrrr@{\hspace{0.4cm}}lrrrrr@{\hspace{0.4cm}}lrrrrr@{\hspace{0.4cm}}lrrrrr}(&1,&1,&1,&1&)\hskip
11.38092pt&(&-1,&1,&1,&1&)\hskip 11.38092pt&(&-1,&-1,&1,&1&)\hskip
11.38092pt&(&-1,&-1,&-1,&1&)\hskip 11.38092pt&(&-1,&-1,&-1,&-1&)\\\
&&&&&\hskip 11.38092pt&(&1,&-1,&1,&1&)\hskip 11.38092pt&(&-1,&1,&-1,&1&)\hskip
11.38092pt&(&-1,&-1,&1,&-1&)\hskip 11.38092pt&&&&&&\\\ &&&&&\hskip
11.38092pt&(&1,&1,&-1,&1&)\hskip 11.38092pt&(&1,&-1,&-1,&1&)\hskip
11.38092pt&(&-1,&1,&-1,&-1&)\hskip 11.38092pt&&&&&&\\\ &&&&&\hskip
11.38092pt&(&1,&1,&1,&-1&)\hskip 11.38092pt&(&-1,&1,&1,&-1&)\hskip
11.38092pt&(&1,&-1,&-1,&-1&)\hskip 11.38092pt&&&&&&\\\ &&&&&\hskip
11.38092pt&&&&&&\hskip 11.38092pt&(&1,&-1,&1,&-1&)\hskip
11.38092pt&&&&&&\hskip 11.38092pt&&&&&&\\\ &&&&&\hskip 11.38092pt&&&&&&\hskip
11.38092pt&(&1,&1,&-1,&-1&)\hskip 11.38092pt&&&&&&\hskip 11.38092pt&&&&&&\\\
\end{array}$

are averaged before performing the ensemble averages. By doing so we are
assuming that the dependence of the form factors on the H4 invariants
$p^{[n]}$ is small or negligible.

##  IV Lattice form factors

Figure 3: Bare gluon propagator in the Landau gauge using the data from the
simulations with the $32^{4}$ and $48^{4}$ lattices. The data reported is a
subset of all the lattice momenta. Indeed, to handle the effects associated
with the breaking of rotational symmetry we follow the recipe of Leinweber
_et al._ (1999) and perform the so-called cylindrical and conical cuts in
momentum space, while for $p<0.7$ GeV all data is included.

Figure 4: Dimensionless bare lattice form factors $F^{(i)}$ for the various
kinematics. Black dots are the data from the simulations on the $32^{4}$
lattice, while the red squares refer to the simulation on $48^{4}$. Due to the
large statistical errors the data from the larger volume is omitted in some of
the plots. Note the different vertical scales used in each Fig. See main text
for details.

The continuum form factors $F^{(i)}$ defined in Eqs (31) to (33) were measured
in the simulations using their lattice implementations as given in App. B. For
the ensembles referred previously, it turns out that only the simulations
using the $32^{4}$ lattice have a good signal-to-noise ratio on a range of
momenta. On the other hand, the results coming from the simulation with the
$48^{4}$ lattice show a good signal-to-noise ratio only for $F^{(0)}$ but not
for the amputated functions. Anyway, despite the relative large errors
associated with this latter lattice, the outcome of the simulation confirms
the tendency observed in the $32^{4}$ results and the two sets of data are
compatible within one standard deviation. In all cases, the statistical errors
were computed with the bootstrap method with a confidence level of 67.5%. The
bare lattice data for the form factors will be described in terms of the
improved lattice momenta

$p_{\mu}=(2/a)\,\sin(\pi\,n_{\mu}/L)\qquad\mbox{ with }\qquad
n_{\mu}=-L/2,\,-L/2+1,\,\cdots,\,0,\,1,\,\cdots,\,L/2-1\ ,$ (43)

where $L$ is the length of the lattice. At least for the gluon propagator the
use of the improved momenta handles part of the lattice effects.

The bare form factors $F^{(0)}$, $F^{(1)}$ and $F^{(2)}$ associated with the
full Green functions, i.e. that incorporate the contributions of the gluon
propagators, are shown in Fig. 4 for all the kinematics and for all the
momenta accessed in the simulation. The black circles are the data from
simulations using the $32^{4}$ lattice volume, while the red squares is the
data measured on the $48^{4}$ lattice volume. Note the different vertical
scales that are associated with each of the $F^{(i)}$. Their order of
magnitude can be understood looking at the bare gluon propagator function
$D(p^{2})$, see Fig. 3, that is included in the data of Fig. 4.

The data reported shows there is a range of momenta where there is a clear
sign for all the form factors that goes up to $p\sim 1$ GeV. For higher
momenta the statistical errors prevent the access to the behaviour of the
various $F^{(i)}$. Reduced errors can be accessed by using larger statistical
ensembles of gauge configurations. However, the data also shows that the
measurement of the four gluon Green functions is possible for gauge ensembles
with, at least, about 10k configurations.

Figure 5: Dimensionless bare lattice amputated form factors $F^{(i)}$ for the
simulation on $32^{4}$ lattice. The data from the $48^{4}$ simulation was
omitted as it has a two large statistical error to add any useful information.
For the lower momenta the data from the simulation with a $48^{4}$, that have
smaller statistical errors, follow the trend observed in the $32^{4}$
simulation and the two sets of data are compatible within one standard
deviation.

Figure 6: Dimensionless bare lattice amputated form factors $F^{(i)}$ as a
function of $s$ for the simulation on $32^{4}$ lattice. The plots include the
data for all the kinematical configurations. In the right plots the data for
kinematical configuration with the large statistical error was omitted for a
better reading on the dependence of the form factors with $s$.

Let us now consider the same form factors but associated with the amputated
Green functions. The form factors describing the 1-PI Green function, i.e. the
amputated Green function, can be computed from the $F^{(i)}$ dividing the
later functions by the gluon propagators that are associated with each of the
external legs of the full Green function. For the size of each gauge ensemble,
the statistical errors associated with the two point function are tiny, see
the results of Fig. 3, and will be ignored when performing the division. Then,
the statistical errors for the amputated Green function form factors come only
from the statistical errors on the $F^{(i)}$.

On a lattice simulation the (bare) computed functions relevant for our purpose
are, using now a simplified notation, $a^{4}\langle A(1)\cdots A(4)\rangle$
for the four point function and $a^{2}\langle A(1)A(2)\rangle$ for the
propagator. From the division of the bare $F^{(i)}$ by the bare the gluon
propagator, the would be amputated Green functions are

$\frac{a^{4}\langle A(1)\cdots A(4)\rangle}{\Big{(}a^{2}\langle
A(1)A(2)\rangle\Big{)}^{4}}=\frac{\langle A(1)\cdots
A(4)\rangle}{a^{4}\,\Big{(}\langle A(1)A(2)\rangle\Big{)}^{4}}\ .$ (44)

To arrive at the corresponding dimensionless amputated Green functions, these
ratios have to be multiplied by $a^{4}$. Then, the bare lattice amputated
Green function reads

$a^{4}~{}\frac{a^{4}\langle A(1)\cdots A(4)\rangle}{\Big{(}a^{2}\langle
A(1)A(2)\rangle\Big{)}^{4}}$ (45)

and can be seen in Fig. 5 as a function of the momenta. The data reported
shows that there is a clear Monte Carlo signal for momenta up to $\sim$ 1 GeV
for all the amputated form factors. For larger momenta the statistical errors
become quite large and it is hard to comment on the behaviour of the amputated
$F^{(i)}$. The dependence on $p$ shows an hierarchy between the various
amputated form factors, with the form factor associated with tree level tensor
structure $F^{(0)}$ being the largest of the three $F^{(i)}$. In general, for
the same form factor, the statistical error associated with kinematical
configuration $(p,\,p,\,p,\,-3p)$ is the largest. This can be understood
looking at the lattice operators given in App. B. Indeed, for this kinematical
configuration the number of averages performed for each configurations is
substantially smaller than for the remaining ones and, therefore, the data
associated with $(p,\,p,\,p,\,-3p)$ should have larger fluctuations, when
compared with the other two kinematics investigated here.

In what concerns the dependence of the amputated form factors with momentum,
the curves suggest that $F^{(0)}$ is essentially constant for $p\gtrsim 0.6$
GeV (although for $p\gtrsim 1$ GeV the large statistical errors prevent any
firm conclusion) and decreases for smaller momenta. The systematics for
$F^{(1)}$ and $F^{(2)}$ are harder to understand from Fig. 5. As discussed
below, by taking into account the Bose symmetry and its implication on the
amputated form factors gives a clear understanding on their dependence with
momentum scales.

The gluon field is a bosonic field and, therefore, the form factors should
depend on momenta as described in Eq. (34). In this sense, the form factors
should be a function not of $p$ but of

$s=\left(p^{2}_{1}+p^{2}_{2}+p^{2}_{3}+p^{2}_{4}\right)/4$ (46)

and of the combination of scalar products of the various momenta mentioned in
Eq. (34). Let us ignore this latter dependence and redo the plots, looking at
the amputated form factors as a function of $s$ and combining the various
kinematics in the same plot. The lattice data for the three amputated form
factors is reported in Fig. 6. In general, the lattice data for the amputated
form factors is compatible with a dependence on $s$. For the level of
statistical precision achieved in the simulation, the description of the
lattice does not seem to require any further variable than $s$ itself. The
data of the amputated $F^{(0)}$ is compatible with a constant for $s\gtrsim 1$
GeV. However, $F^{(0)}$ is slightly suppressed at lower momenta. On the other
hand $F^{(1)}$ and $F^{(2)}$ seem to be constant for $s\gtrsim 1$ GeV but
increase when $s$ becomes smaller. Moreover, the lattice data for the
amputated $F^{(1)}$ suggests that this form factor should change sign at
$s\sim 0.3$ GeV. No change of sign is seen or suggested for $F^{(2)}$. The
relative strong increase in the $F^{(1)}$ and $F^{(2)}$ observed when $s$
approaches zero, seem to suggest that these form factors may have logarithmic
divergences in the IR region. However, given the limited access to the low
momenta region, the question of the logarithmic divergences is not solved
within the current simulation. To answer this question one should have access
to the deep infrared region that require simulations using larger lattices or
with different lattice spacing.

##  V Summary and Conclusions

Herein, the calculation of the four gluon one-particle irreducible Green
function, in the Landau gauge, using lattice QCD simulations is addressed. The
determination of the associated form factors from the full Green function is
discussed and some of the 1-PI Green function form factors are computed. Our
study shows that lattice simulations to access the four gluon 1-PI Green
function require large ensembles of gauge configurations, meaning 10k or more
gauge configurations per ensemble. Moreover, to resolve the 1-PI contribution
to the full four gluon Green function, that is the Green function directly
measured in the simulation, the kinematics of the external legs have to be
carefully chosen. In particular, for the kinematics characterized by a unique
momentum, it is possible to disentangle the four gluon 1-PI contribution from
the disconnected parts and from the diagrams associated with the three gluon
1-PI in the full Green function.

In order to arrive at such large ensembles of configurations, the simulations
reported use either a hypercubic lattice size of $32^{4}$, that has a physical
volume of $(\sim 3.3\mbox{ fm})^{4}$, or a hypercubic lattice size of
$48^{4}$, whose physical volume is about $(\sim 4.9\mbox{ fm})^{4}$. For the
smaller lattice our ensemble has about 10K configurations, for the larger
lattice the size of the ensemble is of the order of 4K. This prevents the
larger lattice volume simulation to achieve a good signal-to-noise ratio and
the statistical errors for the 1-PI form factors for the $48^{4}$ lattice are
quite large. Nevertheless, they confirm the tendency observed in the
simulation with the $32^{4}$ lattice. We are currently working on the
production of larger statistical ensembles to improve our results, i.e. to
achieve a good signal for larger range of momenta and to explore other
kinematical configurations. The results of such simulations will be made
public as soon as possible.

The complete description of the four gluon 1-PI irreducible Green function
calls for a tensor basis with more than one hundred independent operators.
However, for the kinematical configurations under consideration, the number of
relevant operators that contribute to full Green functions is relatively small
and, in some cases, can be handled exactly. Of the possible tensor operators
that define the basis for the 1-PI, we have considered the three operators
(11), (12) and (13), measured the associated form factors $F^{(0)}$,
$F^{(1)}$, $F^{(2)}$ defined in (31), (32), (33), respectively, and their bare
amputated versions that appear in the 1-PI four gluon irreducible Green
function.

Figure 7: Dimensionless bare lattice amputated form factors $F^{(0)}$ as a
function of $s$ for the simulation on $32^{4}$ lattice as in Fig. 6 replacing
the various data points for the same momentum with an weighted average using
the statistical error squared as weight.

The tensor operator used to define the amputated form factor $F^{(0)}$ is the
operator that appears in the tree level Feynman rule for the four gluon
vertex. Our simulation gives an amputated form factor that is essentially
constant for momentum above $\sim 0.5$ GeV (see Fig. 6), with the data
suggesting a suppression at small momenta. Despite the large ensemble, for
momentum above $\sim 1$ GeV it is difficult to disentangle its functional
form. For momenta in the range $\sim 1-1.5$ GeV, despite the size of the
statistical errors, the data seems to be compatible with a constant value;
this can be better seen by replacing the various data points for the same
momentum by their weighted averages using the statistical error squared as a
weight as reported in Fig. 7 for $F^{(0)}$. This result is in qualitative
agreement with the outcome of recent Dyson-Schwinger calculations Binosi  _et
al._ (2014); Cyrol  _et al._ (2015); Colaço  _et al._ (2023) but not with the
calculation described in Kellermann and Fischer (2008).

The analysis of $F^{(1)}$ shows a negative form factor whose absolute value is
about a quarter of the absolute value of $F^{(0)}$, it has the opposite sign
of this latter form factor and seems to approach zero as $s$ is decreased.
This behaviour suggests a change of sign for $F^{(1)}$ at $s\sim 0.15$ GeV2
and suggests a divergent behaviour for $F^{(1)}$ in the deep infrared region.
Again, the observed $F^{(1)}$ is qualitatively in agreement with the Dyson-
Schwinger result reported in Colaço  _et al._ (2023).

The bare amputated $F^{(2)}$ form factor has, apart the possible change of
sign, a similar behaviour as observed for $F^{(1)}$, i.e. $F^{(2)}\sim
F^{(0)}/4$, with both form factors having the same sign and with the $F^{(2)}$
data suggesting a possible $\log$ divergence in the deep infrared region. The
observed lattice data for $F^{(2)}$ is qualitatively in agreement with the
recent Dyson-Schwinger calculation reported in Colaço  _et al._ (2023). Note
that the remaining Dyson-Schwinger studies do not access either $F^{(1)}$ or
$F^{(2)}$.

## Acknowledgments

The authors acknowledge Arlene C. Aguilar, Mauricio N. Ferreira, Joannis
Papavassiliou and Leonardo R. Santos for helpful discussions and careful
reading of the manuscript. This work was partly supported by the FCT –
Fundação para a Ciência e a Tecnologia, I.P., under Projects Nos.
UIDB/04564/2020 (https://doi.org/10.54499/UIDB/04564/2020), UIDP/04564/2020
(https://doi.org/10.54499/UIDP/04564/2020). and CERN/FIS-PAR/0023/2021. P. J.
S. acknowledges financial support from FCT contract CEECIND/00488/2017
(https://doi.org/10.54499/CEECIND/00488/2017/CP1460/CT0030). The authors
acknowledge the Laboratory for Advanced Computing at the University of Coimbra
(http://www.uc.pt/lca) for providing access to the HPC resources that have
contributed to the research within this paper. Access to Navigator was partly
supported by the FCT Advanced Computing Project 2021.09759.CPCA.

##  Appendix A Three gluon contributions for proportional momenta

The three-gluon one-particle diagram can be analysed using the Ball-Chiu
tensor basis introduced in Ball and Chiu (1980); see also Ahmadiniaz and
Schubert (2013). For the momenta configuration of type $p_{1}=p$,
$p_{2}=\eta\,p$ and $p_{3}=-(1+\eta)\,p$, only the so called longitudinal
component contributes. Following the notation of Ball-Chiu, the one-particle
longitudinal component reads

$\displaystyle(1-\eta)\,A(p^{2}_{1},\,p^{2}_{2};\,p^{2}_{3})~{}g_{\mu_{1}\mu_{2}}~{}p_{\mu_{3}}+(1+2\,\eta)\,A(p^{2}_{2},\,p^{2}_{3};\,p^{2}_{1})~{}g_{\mu_{2}\mu_{3}}~{}p_{\mu_{1}}-(2+\,\eta)\,A(p^{2}_{3},\,p^{2}_{1};\,p^{2}_{2})~{}g_{\mu_{3}\mu_{1}}~{}p_{\mu_{2}}$
$\displaystyle\quad+~{}(1+\eta)\,B(p^{2}_{1},\,p^{2}_{2};\,p^{2}_{3})~{}g_{\mu_{1}\mu_{2}}~{}p_{\mu_{3}}-(1+2\,\eta)\,B(p^{2}_{2},\,p^{2}_{3};\,p^{2}_{1})~{}g_{\mu_{2}\mu_{3}}~{}p_{\mu_{1}}-\eta\,B(p^{2}_{3},\,p^{2}_{1};\,p^{2}_{2})~{}g_{\mu_{3}\mu_{1}}~{}p_{\mu_{2}}$
$\displaystyle\quad-~{}\eta(1-\eta)\,C(p^{2}_{1},\,p^{2}_{2};\,p^{2}_{3})~{}p^{2}~{}P^{\perp}_{\mu_{1}\mu_{2}}(p)~{}p_{\mu_{3}}+\eta(1+\,\eta)\,C(p^{2}_{2},\,p^{2}_{3};\,p^{2}_{1})~{}p^{2}~{}P^{\perp}_{\mu_{2}\mu_{3}}(p)~{}p_{\mu_{1}}$
$\displaystyle\quad-(1+\eta)(2+\eta)\,C(p^{2}_{3},\,p^{2}_{1};\,p^{2}_{2})~{}p^{2}~{}P^{\perp}_{\mu_{1}\mu_{3}}(p)~{}p_{\mu_{2}}~{}-~{}2\,\eta\,(1+\eta)S(p^{2}_{1},\,p^{2}_{2},\,p^{2}_{3})~{}p_{\mu_{1}}\,p_{\mu_{2}}\,p_{\mu_{3}}\,$
(47)

where $A$, $B$, $C$, $S$ are the Ball-Chiu form factors and

$P^{\perp}_{\mu\nu}(p)=g_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{p^{2}}\qquad\mbox{
with }\qquad p^{\mu}P^{\perp}_{\mu\nu}(p)=0$ (48)

is the orthogonal projector for vector fields. The external legs are
proportional to $p$ at least for one of the Lorentz indices, given that in the
contribution to the four-gluon connected Green function all the external legs
are contracted with a gluon propagator that, in the Landau gauge, is
proportional to $P^{\perp}_{\mu\nu}(p)$, given that $P^{\perp}_{\mu\nu}(\kappa
p)=P^{\perp}_{\mu\nu}(p)$, where $\kappa$ is a constant, then it follows that
for this class of kinematical configuration the terms in Fig. 2 that contain
the three-gluon irreducible diagram vanish. Note that this result is valid
even if one of the momenta vanishes.

##  Appendix B Implementation of color-Lorentz projectors to measure $F^{(i)}$

The tensor operators that define the lattice form factors $F^{(i)}$, see Eqs
(31) to (33), require the SU(3) totally antisymmetric group structure
constants $f_{abc}$, the totally symmetric structure constants $d_{abc}$ or
color diagonal $\delta^{ab}$ operators. The first two types of operators can
be computed from color traces of the generators given in the fundamental
representation $t^{a}$ as

$\mbox{Tr}\left(t^{a}\,t^{b}\,t^{c}\right)=\frac{1}{4}\big{(}d_{abc}+i\,f_{abc}\big{)}$
(49)

or, alternatively,

$-\,2\,i\,\mbox{Tr}\left(\left[t^{a}\,,\,t^{b}\right]\,t^{c}\right)=f_{abc}\qquad\mbox{
and }\qquad
2\,\,\mbox{Tr}\left(\left\\{t^{a}\,,\,t^{b}\right\\}\,t^{c}\right)=d_{abc}\ .$
(50)

The Euclidean version of the tensors in Eqs (11) to (13) is

$\displaystyle\widetilde{\Gamma}^{(0)}\,{}^{abcd}_{\mu\nu\eta\zeta}$
$\displaystyle=$ $\displaystyle
f_{abr}f_{cdr}\Big{(}\delta_{\mu\eta}\,\delta_{\nu\zeta}-\delta_{\mu\zeta}\,\delta_{\nu\eta}\Big{)}+f_{acr}f_{bdr}\Big{(}\delta_{\mu\nu}\,\delta_{\eta\zeta}-\delta_{\mu\zeta}\,\delta_{\nu\eta}\Big{)}+f_{adr}f_{bcr}\Big{(}\delta_{\mu\nu}\,\delta_{\eta\zeta}-\delta_{\mu\eta}\,\delta_{\nu\zeta}\Big{)}\
,$ (51) $\displaystyle\widetilde{\Gamma}^{(1)}\,{}^{abcd}_{\mu\nu\eta\zeta}$
$\displaystyle=$ $\displaystyle
d_{abr}d_{cdr}\Big{(}\delta_{\mu\eta}\,\delta_{\nu\zeta}+\delta_{\mu\zeta}\,\delta_{\nu\eta}\Big{)}+d_{acr}d_{bdr}\Big{(}\delta_{\mu\zeta}\,\delta_{\nu\eta}+\delta_{\mu\nu}\,\delta_{\eta\zeta}\Big{)}+d_{adr}d_{bcr}\Big{(}\delta_{\mu\nu}\,\delta_{\eta\zeta}+\delta_{\mu\eta}\,\delta_{\nu\zeta}\Big{)}\
,$ (52) $\displaystyle\widetilde{\Gamma}^{(2)}\,{}^{abcd}_{\mu\nu\eta\zeta}$
$\displaystyle=$
$\displaystyle\Big{(}\delta^{ab}\,\delta^{cd}+\delta^{ac}\,\delta^{bd}+\delta^{ad}\,\delta^{bc}\Big{)}~{}\Big{(}\delta_{\mu\nu}\,\delta_{\eta\zeta}+\delta_{\mu\eta}\,\delta_{\nu\zeta}+\delta_{\mu\zeta}\,\delta_{\nu\eta}\Big{)}$
(53)

and

$F^{(i)}~{}=~{}\quad\widetilde{\Gamma}^{(i)}\,{}^{abcd}_{\mu\nu\eta\zeta}\quad\langle
A^{a}_{\mu}(p_{1})~{}A^{b}_{\nu}(p_{2})A^{c}_{\eta}(p_{3})A^{d}_{\zeta}(p_{4})\rangle\
.$ (54)

Then, plugging in these definitions it turns out that for the (K. Conf. 1)
kinematics

$\boxed{p_{1}=0,\,p_{2}=p_{3}=p,~{}p_{4}=-\,2\,p}$ $\displaystyle
F^{(0)}(p^{2})$ $\displaystyle=$
$\displaystyle-\,8~{}\sum_{r}~{}\Bigg{\\{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\mu}(p)\,,\,A_{\nu}(-2p)\right]\,t^{r}\Big{)}$
(55) $\displaystyle\hskip
56.9055pt-\,\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\nu}(p)\,,\,A_{\mu}(-2p)\right]\,t^{r}\Big{)}$
$\displaystyle\hskip
85.35826pt+\,\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(-2p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\mu}(p)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\Bigg{\\}}\
,$ $\displaystyle F^{(1)}(p^{2})$ $\displaystyle=$ $\displaystyle
8~{}\sum_{r}~{}\Bigg{\\{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(-2p)\right\\}\,t^{r}\Big{)}$
(56) $\displaystyle\hskip
85.35826pt+\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\nu}(p)\,,\,A_{\mu}(-2p)\right\\}\,t^{r}\Big{)}$
$\displaystyle\hskip
113.81102pt+\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(-2p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\Bigg{\\}}\
,$ $\displaystyle F^{(2)}(p^{2})$ $\displaystyle=$ $\displaystyle
4~{}\Bigg{\\{}2\,\,\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\mu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\nu}(-2p)\,\Big{)}$
(57) $\displaystyle\hskip
56.9055pt+2\,\,\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(-2p)\,\Big{)}$
$\displaystyle\hskip
85.35826pt+\,2\,\,\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\mu}(-2p)\Big{)}$
$\displaystyle\hskip
113.81102pt+\,2\,\,\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(-2p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(p)\Big{)}$
$\displaystyle\hskip
142.26378pt+\,\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\mu}(-2p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\nu}(p)\Big{)}\Bigg{\\}}\
,$

for the (K. Conf. 2) kinematics

$\boxed{p_{1}=0,~{}p_{2}=p,~{}p_{3}=2\,p,~{}p_{4}=-\,3\,p}$ $\displaystyle
F^{(0)}(p^{2})$ $\displaystyle=$
$\displaystyle\,-\,4\,~{}\sum_{r}~{}\Bigg{\\{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\mu}(2p)\,,\,A_{\nu}(-3p)\right]\,t^{r}\Big{)}$
(58) $\displaystyle\hskip
71.13188pt-~{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\nu}(2p)\,,\,A_{\mu}(-3p)\right]\,t^{r}\Big{)}$
$\displaystyle\hskip
85.35826pt+~{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(2p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\mu}(p)\,,\,A_{\nu}(-3p)\right]\,t^{r}\Big{)}$
$\displaystyle\hskip
99.58464pt-~{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(2p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\nu}(p)\,,\,A_{\mu}(-3p)\right]\,t^{r}\Big{)}$
$\displaystyle\hskip
113.81102pt+~{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(-3p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\mu}(p)\,,\,A_{\nu}(2p)\right]\,t^{r}\Big{)}$
$\displaystyle\hskip
128.0374pt-~{}\mbox{Tr}\Big{(}\left[A_{\mu}(0)\,,\,A_{\nu}(-3p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\nu}(p)\,,\,A_{\mu}(2p)\right]\,t^{r}\Big{)}\Bigg{\\}}\
,$ $\displaystyle F^{(1)}(p^{2})$ $\displaystyle=$
$\displaystyle\,4\,~{}\sum_{r}~{}\Bigg{\\{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(2p)\,,\,A_{\nu}(-3p)\right\\}\,t^{r}\Big{)}$
(59) $\displaystyle\hskip
71.13188pt+~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\nu}(2p)\,,\,A_{\mu}(-3p)\right\\}\,t^{r}\Big{)}$
$\displaystyle\hskip
85.35826pt+~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(2p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(-3p)\right\\}\,t^{r}\Big{)}$
$\displaystyle\hskip
99.58464pt+~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(2p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\nu}(p)\,,\,A_{\mu}(-3p)\right\\}\,t^{r}\Big{)}$
$\displaystyle\hskip
113.81102pt+~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(-3p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(2p)\right\\}\,t^{r}\Big{)}$
$\displaystyle\hskip
128.0374pt+~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(0)\,,\,A_{\nu}(-3p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\nu}(p)\,,\,A_{\mu}(2p)\right\\}\,t^{r}\Big{)}\Bigg{\\}}\
,$ $\displaystyle F^{(2)}(p^{2})$ $\displaystyle=$ $\displaystyle
4~{}\Bigg{\\{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\mu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(2p)\,\,A_{\nu}(-3p)\,\Big{)}$
(60) $\displaystyle\hskip
28.45274pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\mu}(2p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\nu}(-3p)\Big{)}$
$\displaystyle\hskip
42.67912pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\mu}(-3p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\nu}(2p)\Big{)}$
$\displaystyle\hskip
56.9055pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(2p)\,\,A_{\nu}(-3p)\Big{)}$
$\displaystyle\hskip
71.13188pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(2p)\,\,A_{\mu}(-3p)\Big{)}$
$\displaystyle\hskip
85.35826pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(2p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(-3p)\Big{)}$
$\displaystyle\hskip
99.58464pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(2p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\mu}(-3p)\Big{)}$
$\displaystyle\hskip
113.81102pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(-3p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(2p)\Big{)}$
$\displaystyle\hskip
128.0374pt+~{}\mbox{Tr}\Big{(}A_{\mu}(0)\,\,A_{\nu}(-3p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\mu}(2p)\Big{)}\Bigg{\\}}$

and for the (K. Conf. 3) kinematics

$\boxed{p_{1}=p_{2}=p_{3}=p,~{}p_{4}=-\,3\,p}$ $\displaystyle F^{(0)}(p^{2})$
$\displaystyle=$
$\displaystyle\,24~{}\sum_{r}~{}\mbox{Tr}\Big{(}\left[A_{\mu}(p)\,,\,A_{\nu}(p)\right]\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left[A_{\nu}(p)\,,\,A_{\mu}(-3p)\right]\,t^{r}\Big{)}\
,$ (61) $\displaystyle F^{(1)}(p^{2})$ $\displaystyle=$ $\displaystyle
24~{}\sum_{r}~{}\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(p)\right\\}\,t^{r}\Big{)}\,\mbox{Tr}\Big{(}\left\\{A_{\mu}(p)\,,\,A_{\nu}(-3p)\right\\}\,t^{r}\Big{)}\
,$ (62) $\displaystyle F^{(2)}(p^{2})$ $\displaystyle=$ $\displaystyle
12~{}\Bigg{\\{}\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\mu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\nu}(p)\,\,A_{\nu}(-3p)\,\Big{)}$
(63) $\displaystyle\hskip
56.9055pt+\,2\,\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(p)\Big{)}\,\mbox{Tr}\Big{(}A_{\mu}(p)\,\,A_{\nu}(-3p)\Big{)}\Bigg{\\}}\
.$

The expressions in Eqs (55) to (63) will be used in the simulations to measure
each of the lattice form factors $F^{(i)}$.

## References

  * Cucchieri and Mendes (2007) A. Cucchieri and T. Mendes, PoS LATTICE2007, 297 (2007), arXiv:0710.0412 [hep-lat] .
  * Aguilar  _et al._ (2008) A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D 78, 025010 (2008), arXiv:0802.1870 [hep-ph] .
  * Fischer  _et al._ (2009) C. S. Fischer, A. Maas, and J. M. Pawlowski, Annals Phys. 324, 2408 (2009), arXiv:0810.1987 [hep-ph] .
  * Bogolubsky  _et al._ (2009) I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Phys. Lett. B 676, 69 (2009), arXiv:0901.0736 [hep-lat] .
  * Oliveira and Silva (2012) O. Oliveira and P. J. Silva, Phys. Rev. D 86, 114513 (2012), arXiv:1207.3029 [hep-lat] .
  * Ayala  _et al._ (2012) A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, and J. Rodriguez-Quintero, Phys. Rev. D 86, 074512 (2012), arXiv:1208.0795 [hep-ph] .
  * Duarte  _et al._ (2017) A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 96, 098502 (2017), arXiv:1704.02864 [hep-lat] .
  * Alkofer  _et al._ (2004) R. Alkofer, W. Detmold, C. S. Fischer, and P. Maris, Phys. Rev. D 70, 014014 (2004), arXiv:hep-ph/0309077 .
  * Strauss  _et al._ (2012) S. Strauss, C. S. Fischer, and C. Kellermann, Phys. Rev. Lett. 109, 252001 (2012), arXiv:1208.6239 [hep-ph] .
  * Falcão  _et al._ (2020) A. F. Falcão, O. Oliveira, and P. J. Silva, Phys. Rev. D 102, 114518 (2020), arXiv:2008.02614 [hep-lat] .
  * Falcão and Oliveira (2022) A. F. Falcão and O. Oliveira, Phys. Rev. D 106, 114022 (2022), arXiv:2209.14815 [hep-lat] .
  * Boito  _et al._ (2023) D. Boito, A. Cucchieri, C. Y. London, and T. Mendes, JHEP 02, 144 (2023), arXiv:2210.10490 [hep-lat] .
  * Ilgenfritz  _et al._ (2007) E. M. Ilgenfritz, M. Muller-Preussker, A. Sternbeck, A. Schiller, and I. L. Bogolubsky, Braz. J. Phys. 37, 193 (2007), arXiv:hep-lat/0609043 .
  * Schleifenbaum  _et al._ (2005) W. Schleifenbaum, A. Maas, J. Wambach, and R. Alkofer, Phys. Rev. D 72, 014017 (2005), arXiv:hep-ph/0411052 .
  * Cucchieri  _et al._ (2008) A. Cucchieri, A. Maas, and T. Mendes, Phys. Rev. D 77, 094510 (2008), arXiv:0803.1798 [hep-lat] .
  * Huber and von Smekal (2013) M. Q. Huber and L. von Smekal, JHEP 04, 149 (2013), arXiv:1211.6092 [hep-th] .
  * Binosi  _et al._ (2013) D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 87, 125026 (2013), arXiv:1304.2594 [hep-ph] .
  * Boucaud  _et al._ (2014) P. Boucaud, M. Brinet, F. De Soto, V. Morenas, O. Pène, K. Petrov, and J. Rodríguez-Quintero, JHEP 04, 086 (2014), arXiv:1310.4087 [hep-ph] .
  * Eichmann  _et al._ (2014) G. Eichmann, R. Williams, R. Alkofer, and M. Vujinovic, Phys. Rev. D 89, 105014 (2014), arXiv:1402.1365 [hep-ph] .
  * Athenodorou  _et al._ (2016) A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, Phys. Lett. B 761, 444 (2016), arXiv:1607.01278 [hep-ph] .
  * Duarte  _et al._ (2016a) A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 074502 (2016a), arXiv:1607.03831 [hep-lat] .
  * Rodríguez-Quintero  _et al._ (2017) J. Rodríguez-Quintero, A. Athenodorou, D. Binosi, P. Boucaud, F. de Soto, J. Papavassiliou, and S. Zafeiropoulos, EPJ Web Conf. 137, 03018 (2017).
  * Binosi and Papavassiliou (2018) D. Binosi and J. Papavassiliou, Phys. Rev. D 97, 054029 (2018), arXiv:1709.09964 [hep-ph] .
  * Vujinovic and Mendes (2019) M. Vujinovic and T. Mendes, Phys. Rev. D 99, 034501 (2019), arXiv:1807.03673 [hep-lat] .
  * Aguilar  _et al._ (2019) A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D 99, 094010 (2019), arXiv:1903.01184 [hep-ph] .
  * Aguilar  _et al._ (2020) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, J. Rodríguez-Quintero, and S. Zafeiropoulos, Eur. Phys. J. C 80, 154 (2020), arXiv:1912.12086 [hep-ph] .
  * Maas and Vujinović (2022) A. Maas and M. Vujinović, SciPost Phys. Core 5, 019 (2022), arXiv:2006.08248 [hep-lat] .
  * Aguilar  _et al._ (2021) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, and J. Rodríguez-Quintero, Phys. Lett. B 818, 136352 (2021), arXiv:2102.04959 [hep-ph] .
  * Figueroa and Peláez (2022) F. Figueroa and M. Peláez, Phys. Rev. D 105, 094005 (2022), arXiv:2110.09561 [hep-th] .
  * Catumba  _et al._ (2022) G. T. R. Catumba, O. Oliveira, and P. J. Silva, EPJ Web Conf. 258, 02008 (2022), arXiv:2111.10312 [hep-lat] .
  * Barrios  _et al._ (2022) N. Barrios, M. Peláez, and U. Reinosa, Phys. Rev. D 106, 114039 (2022), arXiv:2207.10704 [hep-ph] .
  * Aguilar  _et al._ (2023a) A. C. Aguilar, M. N. Ferreira, J. Papavassiliou, and L. R. Santos, Eur. Phys. J. C 83, 549 (2023a), arXiv:2305.05704 [hep-ph] .
  * Pinto-Gómez  _et al._ (2023) F. Pinto-Gómez, F. De Soto, M. N. Ferreira, J. Papavassiliou, and J. Rodríguez-Quintero, Phys. Lett. B 838, 137737 (2023), arXiv:2208.01020 [hep-ph] .
  * Pinto-Gómez (2023) F. Pinto-Gómez, Few Body Syst. 64, 43 (2023).
  * Aguilar  _et al._ (2023b) A. C. Aguilar, M. N. Ferreira, B. M. Oliveira, J. Papavassiliou, and L. R. Santos, Eur. Phys. J. C 83, 889 (2023b), arXiv:2306.16283 [hep-ph] .
  * Davydychev  _et al._ (2001) A. I. Davydychev, P. Osland, and L. Saks, Phys. Rev. D 63, 014022 (2001), arXiv:hep-ph/0008171 .
  * Bender  _et al._ (2002) A. Bender, W. Detmold, C. D. Roberts, and A. W. Thomas, Phys. Rev. C 65, 065203 (2002), arXiv:nucl-th/0202082 .
  * Skullerud and Kizilersu (2002) J. Skullerud and A. Kizilersu, JHEP 09, 013 (2002), arXiv:hep-ph/0205318 .
  * Skullerud  _et al._ (2003) J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, JHEP 04, 047 (2003), arXiv:hep-ph/0303176 .
  * Bhagwat and Tandy (2004) M. S. Bhagwat and P. C. Tandy, Phys. Rev. D 70, 094039 (2004), arXiv:hep-ph/0407163 .
  * Kizilersu  _et al._ (2007) A. Kizilersu, D. B. Leinweber, J.-I. Skullerud, and A. G. Williams, Eur. Phys. J. C 50, 871 (2007), arXiv:hep-lat/0610078 .
  * Matevosyan  _et al._ (2007) H. H. Matevosyan, A. W. Thomas, and P. C. Tandy, J. Phys. G 34, 2153 (2007), arXiv:0706.2393 [nucl-th] .
  * Alkofer  _et al._ (2009) R. Alkofer, C. S. Fischer, F. J. Llanes-Estrada, and K. Schwenzer, Annals Phys. 324, 106 (2009), arXiv:0804.3042 [hep-ph] .
  * He (2009) H.-x. He, Phys. Rev. D 80, 016004 (2009), arXiv:0906.2834 [hep-ph] .
  * Rojas  _et al._ (2013) E. Rojas, J. P. B. C. de Melo, B. El-Bennich, O. Oliveira, and T. Frederico, JHEP 10, 193 (2013), arXiv:1306.3022 [hep-ph] .
  * Aguilar  _et al._ (2014) A. C. Aguilar, D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 90, 065027 (2014), arXiv:1405.3506 [hep-ph] .
  * Peláez  _et al._ (2015) M. Peláez, M. Tissier, and N. Wschebor, Phys. Rev. D 92, 045012 (2015), arXiv:1504.05157 [hep-th] .
  * Binosi  _et al._ (2017) D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 95, 031501 (2017), arXiv:1609.02568 [nucl-th] .
  * Bermudez  _et al._ (2017) R. Bermudez, L. Albino, L. X. Gutiérrez-Guerrero, M. E. Tejeda-Yeomans, and A. Bashir, Phys. Rev. D 95, 034041 (2017), arXiv:1702.04437 [hep-ph] .
  * Oliveira  _et al._ (2018) O. Oliveira, T. Frederico, W. de Paula, and J. P. B. C. de Melo, Eur. Phys. J. C 78, 553 (2018), arXiv:1807.00675 [hep-ph] .
  * Oliveira  _et al._ (2019) O. Oliveira, W. de Paula, T. Frederico, and J. P. B. C. de Melo, Eur. Phys. J. C 79, 116 (2019), arXiv:1807.10348 [hep-ph] .
  * Sultan  _et al._ (2021) M. A. Sultan, K. Raya, F. Akram, A. Bashir, and B. Masud, Phys. Rev. D 103, 054036 (2021), arXiv:1810.01396 [nucl-th] .
  * Kızılersü  _et al._ (2021) A. Kızılersü, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Phys. Rev. D 103, 114515 (2021), arXiv:2103.02945 [hep-lat] .
  * Gao  _et al._ (2021) F. Gao, J. Papavassiliou, and J. M. Pawlowski, Phys. Rev. D 103, 094013 (2021), arXiv:2102.13053 [hep-ph] .
  * Skullerud  _et al._ (2022) J.-I. Skullerud, A. Kızılersü, O. Oliveira, P. Silva, and A. Sternbeck, PoS LATTICE2021, 305 (2022), arXiv:2111.13455 [hep-lat] .
  * El-Bennich  _et al._ (2022) B. El-Bennich, F. E. Serna, R. C. da Silveira, L. A. F. Rangel, A. Bashir, and E. Rojas, Rev. Mex. Fis. Suppl. 3, 0308092 (2022), arXiv:2201.04144 [hep-ph] .
  * Marques  _et al._ (2023) J. Marques, G. Kalusche, T. Mendes, P. J. Silva, J.-I. Skullerud, and O. Oliveira, PoS LATTICE2022, 280 (2023), arXiv:2301.10607 [hep-lat] .
  * Aguilar  _et al._ (2023c) A. C. Aguilar, M. . N. Ferreira, D. Ibañez, and J. Papavassiliou, Eur. Phys. J. C 83, 967 (2023c), arXiv:2308.16297 [hep-ph] .
  * Eichmann  _et al._ (2021) G. Eichmann, J. M. Pawlowski, and J. a. M. Silva, Phys. Rev. D 104, 114016 (2021), arXiv:2107.05352 [hep-ph] .
  * Kellermann and Fischer (2008) C. Kellermann and C. S. Fischer, Phys.Rev. D78, 025015 (2008), arXiv:0801.2697 [hep-ph] .
  * Binosi  _et al._ (2014) D. Binosi, D. Ibañez, and J. Papavassiliou, JHEP 09, 059 (2014), arXiv:1407.3677 [hep-ph] .
  * Cyrol  _et al._ (2015) A. K. Cyrol, M. Q. Huber, and L. von Smekal, Eur. Phys. J. C 75, 102 (2015), arXiv:1408.5409 [hep-ph] .
  * Gracey (2014) J. A. Gracey, Phys. Rev. D 90, 025011 (2014), arXiv:1406.1618 [hep-ph] .
  * Eichmann  _et al._ (2015) G. Eichmann, C. S. Fischer, and W. Heupel, Phys. Rev. D 92, 056006 (2015), arXiv:1505.06336 [hep-ph] .
  * Colaço  _et al._ (2023) M. Colaço, O. Oliveira, and P. J. Silva, in  _40th International Symposium on Lattice Field Theory_ (2023) arXiv:2311.17188 [hep-lat] .
  * Catumba (2021) G. T. R. Catumba, (2021), arXiv:2101.06074 [hep-lat] .
  * Boucaud  _et al._ (2017) P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D 96, 098501 (2017), arXiv:1704.02053 [hep-lat] .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Nucl. Phys. B Proc. Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 .
  * Pippig (2013) M. Pippig, SIAM J. Sci. Comput. 35, C213 (2013).
  * Duarte  _et al._ (2016b) A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 014502 (2016b), arXiv:1605.00594 [hep-lat] .
  * Dudal  _et al._ (2018) D. Dudal, O. Oliveira, and P. J. Silva, Annals Phys. 397, 351 (2018), arXiv:1803.02281 [hep-lat] .
  * Becirevic  _et al._ (2000) D. Becirevic, P. Boucaud, J. P. Leroy, J. Micheli, O. Pene, J. Rodriguez-Quintero, and C. Roiesnel, Phys. Rev. D 61, 114508 (2000), arXiv:hep-ph/9910204 .
  * de Soto and Roiesnel (2007) F. de Soto and C. Roiesnel, JHEP 09, 007 (2007), arXiv:0705.3523 [hep-lat] .
  * Catumba  _et al._ (2021) G. T. R. Catumba, O. Oliveira, and P. J. Silva, Phys. Rev. D 103, 074501 (2021), arXiv:2101.04978 [hep-lat] .
  * Leinweber  _et al._ (1999) D. B. Leinweber, J. I. Skullerud, A. G. Williams, and C. Parrinello (UKQCD), Phys. Rev. D 60, 094507 (1999), [Erratum: Phys.Rev.D 61, 079901 (2000)], arXiv:hep-lat/9811027 .
  * Ball and Chiu (1980) J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2550 (1980), [Erratum: Phys.Rev.D 23, 3085 (1981)].
  * Ahmadiniaz and Schubert (2013) N. Ahmadiniaz and C. Schubert, Nucl. Phys. B 869, 417 (2013), arXiv:1210.2331 [hep-ph] .