#  Schwinger displacement of the quark-gluon vertex

A. C. Aguilar  University of Campinas - UNICAMP, Institute of Physics Gleb
Wataghin,
13083-859 Campinas, São Paulo, Brazil  M. N. Ferreira  Department of
Theoretical Physics and IFIC, University of Valencia and CSIC, E-46100,
Valencia, Spain  D. Ibañez  University Centre EDEM, Muelle de la Aduana, La
Marina de Valencia, 46024, Valencia, Spain.  J. Papavassiliou  Department of
Theoretical Physics and IFIC, University of Valencia and CSIC, E-46100,
Valencia, Spain

###### Abstract

The action of the Schwinger mechanism in pure Yang-Mills theories endows
gluons with an effective mass, and, at the same time, induces a measurable
displacement to the Ward identity satisfied by the three-gluon vertex. In the
present work we turn to Quantum Chromodynamics with two light quark flavors,
and explore the appearance of this characteristic displacement at the level of
the quark-gluon vertex. When the Schwinger mechanism is activated, this vertex
acquires massless poles, whose momentum-dependent residues are determined by a
set of coupled integral equations. The main effect of these residues is to
displace the Ward identity obeyed by the pole-free part of the vertex, causing
modifications to its form factors, and especially the one associated with the
tree-level tensor. The comparison between the available lattice data for this
form factor and the Ward identity prediction reveals a marked deviation, which
is completely compatible with the theoretical expectation for the attendant
residue. This analysis corroborates further the self-consistency of this mass-
generating scenario in the general context of real-world strong interactions.

##  I Introduction

In recent years, the generation of a mass scale in the gauge sector of Quantum
Chromodynamics (QCD) Marciano and Pagels (1978) through the action of the
Schwinger mechanism (SM) Schwinger (1962a, b) has emerged as a particularly
appealing and robust notion Smit (1974); Eichten and Feinberg (1974); Cornwall
(1982); Aguilar  _et al._ (2012a); Binosi  _et al._ (2012); Aguilar  _et al._
(2016a); Eichmann  _et al._ (2021); Aguilar  _et al._ (2022, 2023a);
Papavassiliou (2022); Ferreira and Papavassiliou (2023); for alternative
approaches see Fischer  _et al._ (2009); Wilson  _et al._ (1994); Kondo
(2001); Epple  _et al._ (2008); Campagnari and Reinhardt (2010); Tissier and
Wschebor (2010); Alkofer and Alkofer (2011); Serreau and Tissier (2012);
Siringo (2016); Cyrol  _et al._ (2016); Glazek  _et al._ (2017); Huber
(2020a); Horak  _et al._ (2022). In its most general formulation, the SM is
based on the key observation that, even if a mass is forbidden at the level of
the fundamental Lagrangian, a gauge boson may become massive if its vacuum
polarization develops a pole at zero momentum transfer (“massless pole”). Such
massless poles may occur for a variety of reasons, depending on the
dimensionality of space-time and the dynamical details of each theory Zumino
(1965); Jackiw and Johnson (1973); Jackiw (1973); Cornwall and Norton (1973).

The non-Abelian version of the SM operating in Yang-Mills theories proceeds
via the inclusion of massless scalar excitations in the fundamental vertices
of the theory Smit (1974); Eichten and Feinberg (1974); Poggio  _et al._
(1975); Cornwall (1979, 1982); Aguilar  _et al._ (2008, 2012a), thus altering
profoundly their analytic properties. These excitations are longitudinally
coupled, and arise dynamically, as color-carrying bound states of two gluons
or a ghost-antighost pair Aguilar  _et al._ (2012a); Papavassiliou (2022);
Binosi and Papavassiliou (2018); Aguilar  _et al._ (2023a); Ferreira and
Papavassiliou (2023); Aguilar  _et al._ (2022, 2018); Eichmann  _et al._
(2021). Thanks to the functional equations that couple propagators and
vertices Roberts and Williams (1994); Cloet and Roberts (2014); Alkofer and
von Smekal (2001); Fischer (2006); Roberts (2008); Binosi and Papavassiliou
(2009); Aguilar  _et al._ (2016b); Huber (2020b); Aguilar  _et al._ (2023a);
Papavassiliou (2022); Ferreira and Papavassiliou (2023); Fischer  _et al._
(2009); Eichmann  _et al._ (2021); Aguilar  _et al._ (2022), these structures
make their way to the gluon self-energy, providing the kinematic basis for the
appearance of a mass Aguilar  _et al._ (2012a, 2018); Binosi and Papavassiliou
(2018).

Especially important in this context is the residue function associated with
the pole of the three-gluon vertex, denoted by ${\mathbb{C}}(r^{2})$ Aguilar
_et al._ (2022, 2023a); Papavassiliou (2022); Ferreira and Papavassiliou
(2023). This function acts as the Bethe-Salpeter (BS) amplitude that controls
the formation of the aforementioned bound states out of a pair of gluons
Aguilar  _et al._ (2012a, 2018); Binosi and Papavassiliou (2018). In addition,
it is intimately connected with a characteristic displacement of the Ward
identity (WI) satisfied by the pole-free part of the three-gluon vertex,
namely the component measured on the lattice Athenodorou  _et al._ (2016);
Duarte  _et al._ (2016); Boucaud  _et al._ (2017); Pinto-Gómez  _et al._
(2023); Aguilar  _et al._ (2020, 2021a, 2021b). In particular, the form factor
of the three-gluon vertex associated with the so-called “soft-gluon limit” is
shown to deviate from the WI prediction – the correct result when the SM is
inactive – precisely by an amount ${\mathbb{C}}(r^{2})$. Most importantly, as
was established in Aguilar  _et al._ (2023a, 2022), this displacement is
clearly observed in the lattice data: the ${\mathbb{C}}(r^{2})$ reconstructed
is unequivocally nonvanishing, and in excellent agreement with the BS
prediction Aguilar  _et al._ (2022).

In the present work we consider for the first time the appearance of an
analogous effect at the level of the fully-dressed quark-gluon vertex, thus
implementing the pending generalization of the SM from quarkless QCD (pure
Yang-Mills theory) to the case of QCD with two light quark flavors
($N_{f}=2$). We emphasize that the quark-gluon vertex is a crucial component
of the QCD dynamics, and has been studied extensively from the perturbative
point of view Ball and Chiu (1980); Kizilersu  _et al._ (1995); Davydychev
_et al._ (2001); Gracey (2014, 2011); Bermudez  _et al._ (2017), by means of
continuous non-perturbative methods Bender  _et al._ (1996); Bhagwat and Tandy
(2004); Llanes-Estrada  _et al._ (2006); Holl  _et al._ (2005); Matevosyan
_et al._ (2007); Fischer (2006); Aguilar and Papavassiliou (2011); Qin  _et
al._ (2013); Binosi  _et al._ (2017a); Hopfer  _et al._ (2012); Aguilar  _et
al._ (2012b); Rojas  _et al._ (2013); Williams (2015); Williams  _et al._
(2016); Sanchis-Alepuz and Williams (2015); Peláez  _et al._ (2015); Alkofer
_et al._ (2009); Mitter  _et al._ (2015); Cyrol  _et al._ (2018); Aguilar  _et
al._ (2014); Gao  _et al._ (2021); Binosi  _et al._ (2017a); Aguilar  _et al._
(2017); Oliveira  _et al._ (2019a); Albino  _et al._ (2019); Tang  _et al._
(2019); Albino  _et al._ (2021), and through numerous lattice simulations
Skullerud  _et al._ (2002); Skullerud and Kizilersu (2002); Skullerud  _et
al._ (2003, 2005); Lin (2006); Kızılersü  _et al._ (2021); Kizilersu  _et al._
(2007); Sternbeck  _et al._ (2017); Skullerud  _et al._ (2022); Oliveira  _et
al._ (2016, 2018).

The central result of the present work may be summarized by stating that the
action of the SM endows the quark-gluon vertex with a nontrivial pole content.
The general pole structure admitted by Lorentz invariance and charge
conjugation symmetry is comprised by three form factors, which, in the soft-
gluon limit, are the quark-gluon vertex counterparts of the
${\mathbb{C}}(r^{2})$. In this limit, the essential dynamics are described by
a set of coupled BS equations, which account for the fact that the inclusion
of active quarks permits the additional formation of composite scalars out of
quark-antiquark pairs.

The numerical treatment of these equations demonstrates the nonvanishing
nature of the residue functions associated with the poles of the quark-gluon
vertex. Moreover, the modifications produced by their presence to the
unquenched version of ${\mathbb{C}}(r^{2})$ are rather minimal, and the mass-
generating aspects of the SM remain practically unaltered. These findings
clearly validate the expectation Aguilar  _et al._ (2023a) that the gluon mass
generation in the presence of dynamical quarks is qualitatively similar to the
pure Yang-Mills case. In fact, our preliminary analysis is compatible with the
lattice results of Ayala  _et al._ (2012); Binosi  _et al._ (2017b), which
indicate that the unquenched gluon mass is larger than the quenched one.

Interestingly enough, however, the most important effect of these novel
residue functions is the distinctive displacement that they produce to the WI
satisfied by the pole-free part of the quark-gluon vertex. Focusing on the
form factor associated with the classical tensorial structure,
$\gamma_{\alpha}$, this displacement introduces a sharp difference between the
result computed on the lattice, $\lambda_{1}(p^{2})$, and the prediction of
the WI in the absence of the SM, $\lambda_{1}^{\star}(p^{2})$.

The above predictions may be tested by appealing to the available lattice data
for this form factor Kızılersü  _et al._ (2021). In particular, we focus on
the two asymmetric lattices, denominated “L08” and “L07” in Kızılersü  _et
al._ (2021), which have superior statistics and small current quark masses.
The detailed comparison between $\lambda_{1}(p^{2})$ and
$\lambda_{1}^{\star}(p^{2})$ produces a clear signal of $8\sigma$ and
$6\sigma$, respectively, as shown in Fig. 10. In fact, the resulting curves
for the displacement function are in very good agreement with the BS
prediction. We consider these novel results as an additional indication of the
operation of the SM in QCD.

The article is organized as follows. In Sec. II we review the most salient
features of the SM in the context of a Yang-Mills theory, focusing on the
displacement function associated with the three-gluon vertex, which serves as
the prototype for the ensuing construction. In Sec. III we demonstrate how the
action of the SM modifies the WI prediction for the form factors of the quark-
gluon vertex in the soft-gluon limit, and introduce the corresponding
displacement functions. In Sec. IV we solve the system of coupled BS equations
that determines the set of displacement functions associated with the three-
gluon and quark-gluon vertices. Then, in Sec. V we present the central result
of the present work. In particular, we determine the displacement function
associated with the tree-level component of the quark-gluon vertex from the
difference between the lattice data and the WI prediction. In Sec. VI we
present our discussion and conclusions. Finally, in three Appendices we
provide plots and numerical fits for the lattice inputs used in our
computations, explain the implementation of the necessary transitions between
the three main renormalization schemes employed in the literature, and give
some technical details on the Schwinger-Dyson equation (SDE) determination of
a special form factor.

##  II Schwinger mechanism and Ward identity displacement

In this section we review the main concepts of the WI displacement of the
three-gluon vertex, which is induced by the action of the SM. These
considerations set the stage for accomplishing the main objective of this
work, namely study of the same effect in the case of the quark-gluon vertex.

In the Landau gauge that we use throughout this work, the gluon propagator,
$\Delta^{ab}_{\mu\nu}(q)$, has the form111We work in Minkowski space, passing
key results to Euclidean space for the purposes of numerical evaluation or
comparison with the lattice.
$\Delta^{ab}_{\mu\nu}(q)=-i\delta^{ab}\Delta_{\mu\nu}(q)$, with

$\Delta_{\mu\nu}(q)=\Delta(q^{2}){P}_{\mu\nu}(q)\,,\qquad{P}_{\mu\nu}(q):=g_{\mu\nu}-q_{\mu}q_{\nu}/{q^{2}}\,.$
(1)

The momentum evolution of $\Delta(q^{2})$ in the continuum is determined by
the corresponding SDE, given by

$\Delta^{-1}(q^{2})P_{\mu\nu}(q)=q^{2}P_{\mu\nu}(q)+i\,{\Pi}_{\mu\nu}(q)\,,$
(2)

where $\Pi_{\mu\nu}(q)=q^{2}{\bf\Pi}(q^{2})P_{\mu\nu}(q)$ is the transverse
gluon self-energy, shown diagrammatically in Fig. 1.

Figure 1:  Diagrammatic representation of the fully-dressed unquenched gluon
self-energy, $\Pi_{\mu\nu}(q)$, entering in the definition of the SDE for the
gluon propagator, given by Eq. (2). White (colored) circles indicate fully-
dressed propagators (vertices).

A pivotal component of this SDE is the fully-dressed three-gluon vertex,
${\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,k)$, which may be cast in the form

${\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,k)=gf^{abc}{\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)\,,$
(3)

where $g$ is the gauge coupling, $f^{abc}$ are the SU(3) structure constants,
and $q+r+k=0$.

It turns out that the infrared saturation of the gluon propagator hinges
crucially on the special analytic structure of
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}$. To appreciate this point, recall that,
according to the SM, if ${\bf\Pi}(q^{2})$ acquires a pole with a positive
residue (Euclidean space), $\lim_{q^{2}\to 0}{\bf\Pi}(q^{2})=m^{2}/q^{2}$,
then $\Delta^{-1}(0)=m^{2}$, i.e., the gluon propagator picks up a mass
Schwinger (1962a, b); Eichten and Feinberg (1974); Smit (1974); Poggio  _et
al._ (1975); Cornwall (1979, 1982); Jackiw and Johnson (1973); Jackiw (1973);
Cornwall and Norton (1973); Aguilar  _et al._ (2008, 2022, 2023a); Ferreira
and Papavassiliou (2023); Papavassiliou (2022).
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}$ provides precisely this pole222The ghost-
gluon and four-gluon vertices also develop poles, but their influence on the
issues studied here is rather limited Aguilar  _et al._ (2018, 2022)., as a
result of specific nonperturbative dynamics: composite scalar excitations are
formed (e.g., out of a pair of gluons), which carry color, are longitudinally
coupled, and have vanishing mass Eichten and Feinberg (1974); Smit (1974);
Poggio  _et al._ (1975); Cornwall (1979, 1982); Aguilar  _et al._ (2008,
2012a); Papavassiliou (2022); Aguilar  _et al._ (2022, 2023a). As a result,
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}$ may be cast in the form (see upper row of
Fig. 2)

${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)=\Gamma_{\alpha\mu\nu}(q,r,k)+V_{\alpha\mu\nu}(q,r,k)\,,$
(4)

where $\Gamma_{\alpha\mu\nu}(q,r,k)$ denotes the pole-free part, while
$V_{\alpha\mu\nu}(q,r,k)$ contains longitudinal poles of the type
$q^{\alpha}/q^{2}$, $r^{\mu}/r^{2}$, $k^{\nu}/k^{2}$, and products thereof
Aguilar  _et al._ (2023b). The longitudinal nature of
$V_{\alpha\mu\nu}(q,r,k)$ is a direct consequence of the general form of the
gluon-scalar transition amplitude $I_{\alpha}(q)$ [see Fig. 2], namely
$I_{\alpha}(q)=q_{\alpha}I(q^{2})$, imposed by Lorentz symmetry Aguilar  _et
al._ (2012a); it leads to the important property
${P}_{\alpha^{\prime}}^{\alpha}(q){P}_{\mu^{\prime}}^{\mu}(r){P}_{\nu^{\prime}}^{\nu}(k)V_{\alpha\mu\nu}(q,r,k)=0$.
This last relation eliminates the poles from lattice simulations of
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}$ in the Landau gauge; nonetheless,
characteristic finite effects survive, which manifest themselves as
displacements of certain key quantities.

Focusing on the leading pole in the $q$-channel, we have

$V_{\alpha\mu\nu}(q,r,k)=\frac{q_{\alpha}}{q^{2}}g_{\mu\nu}V_{1}(q,r,k)+\cdots\,,$
(5)

where the ellipsis absorbs all remaining pole contributions, which get
annihilated upon the contraction by the projectors $P^{\mu^{\prime}\mu}(r)$ or
$P^{\nu^{\prime}\nu}(k)$.

Figure 2:  Diagrammatic representation of Eq. (4) (top row), and of Eq. (17)
(bottom row). The black circle indicates the gluon-scalar transition function.

The kinematic limit relevant in this study is that of $q\to 0$. A detailed
analysis Aguilar  _et al._ (2016a) shows that $V_{1}(0,r,-r)=0$, and the
Taylor expansion of $V_{1}(q,r,k)$ around $q=0$ yields

$\lim_{q\to 0}V_{1}(q,r,k)=2(q\cdot r)\,{\mathbb{C}}(r^{2})+{\cal
O}(q^{2})\,,$ (6)

where ${\mathbb{C}}(r^{2}):=[\partial V_{1}(q,r,k)/\partial k^{2}]_{q=0}$.

The function ${\mathbb{C}}(r^{2})$ is of central importance in this approach,
playing three key roles:

(i) ${\mathbb{C}}(r^{2})$ is the BS amplitude describing the bound-state
formation of a massless colored scalar out of a pair of gluons. In particular,
${\mathbb{C}}(r^{2})$ is obtained as the nontrivial solution of an appropriate
linear and homogeneous BS equation, deduced from the SDE of
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}$ in the limit $q\to 0$.

(ii) The pole $1/q^{2}$ appearing in Eq. (5) is transmitted to
${\bf\Pi}(q^{2})$ via the gluon SDE of Eq. (2), triggering the SM. The gluon
mass is then determined by integrals involving ${\mathbb{C}}(r^{2})$, of the
general form333 In terms of the gluon-scalar transition amplitude given by the
exact relation $m^{2}=g^{2}I^{2}(0)$ Smit (1974); Eichten and Feinberg (1974);
Ibañez and Papavassiliou (2013), where $I(0)$ is expressed in terms of
integrals such as the one on the r.h.s. of Eq. (7).

$m^{2}\sim\int\\!\\!d^{4}k\,k^{2}\Delta^{2}(k^{2})\,{\mathbb{C}}(k^{2})\,.$
(7)

(iii) If we define the projector ${\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k):=P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(k)$,
the Slavnov-Taylor identity (STI) Taylor (1971); Slavnov (1972) satisfied by
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)$ may be cast in the form

$q^{\alpha}{\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)\,{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k)\\!=\\!F(q^{2})\,[\Delta^{-1}(k^{2})H_{\nu\mu}(k,q,r)\\!-\\!\Delta^{-1}(r^{2})H_{\mu\nu}(r,q,k)]\,{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k)\,,$ (8)

where $F(q^{2})$ is the ghost dressing function, and $H_{\nu\mu}(k,q,r)$ the
ghost-gluon kernel. Note that, by virtue of Eqs. (4) and (5), the l.h.s. of
Eq. (8) becomes

$q^{\alpha}{\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)\,{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k)=[q^{\alpha}\Gamma_{\alpha\mu\nu}(q,r,k)+g_{\mu\nu}V_{1}(q,r,k)]{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k)\,.$ (9)

When the limit $q\to 0$ of both sides of Eq. (8) is taken, one obtains the
corresponding WI, which is tantamount to the so-called “soft-gluon limit” of
the three-gluon vertex. In doing so, Eq. (6) is triggered and
${\mathbb{C}}(r^{2})$ makes its appearance. Moreover, as has been shown in
detail in Aguilar  _et al._ (2021b),

${\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,-r)\Gamma_{\alpha\mu\nu}(0,r,-r)=2r_{\alpha}P_{\mu^{\prime}\nu^{\prime}}(r)\mathit{L}_{{sg}}(r^{2})\,,$
(10)

where $\mathit{L}_{{sg}}(r^{2})$ is precisely the form factor measured on the
lattice through the ratio

$\mathit{L}_{{sg}}(r)=\left.\frac{\Gamma_{0}^{\alpha\mu\nu}(q,r,k)P_{\alpha\alpha^{\prime}}(q)P_{\mu\mu^{\prime}}(r)P_{\nu\nu^{\prime}}(k){\rm{I}}\\!\Gamma^{\alpha^{\prime}\mu^{\prime}\nu^{\prime}}(q,r,k)}{\Gamma_{0}^{\alpha\mu\nu}(q,r,k)P_{\alpha\alpha^{\prime}}(q)P_{\mu\mu^{\prime}}(r)P_{\nu\nu^{\prime}}(k)\Gamma_{0}^{\alpha^{\prime}\mu^{\prime}\nu^{\prime}}(q,r,k)}\right|_{q\to
0}\,,$ (11)

where $\Gamma_{0}^{\alpha\mu\nu}$ is the tree-level expression of the three-
gluon vertex, given by

$\Gamma_{\\!0}^{\alpha\mu\nu}(q,r,k)=(q-r)^{\nu}g^{\alpha\mu}+(r-k)^{\alpha}g^{\mu\nu}+(k-q)^{\mu}g^{\alpha\nu}\,.$
(12)

The final upshot of these considerations is that the function
${\mathbb{C}}(r^{2})$ leads to the displacement of the soft-gluon limit
according to (Euclidean space)

$\underbrace{\mathit{L}_{{sg}}(r^{2})}_{{\rm
SM\,on}}\,=\,\underbrace{\mathit{L}_{{sg}}^{\\!\star}(r^{2})}_{\rm
SM\,off}\,\,+\underbrace{{\mathbb{C}}(r^{2})}_{\begin{subarray}{c}{\rm
displacement}\\\ {\rm function}\end{subarray}}\,,$ (13)

where

$\mathit{L}_{{sg}}^{\\!\star}(r^{2})=[{\rm
r.h.s.\,\,of\,\,STI\,\,in\,\,Eq.\leavevmode\nobreak\ (\ref{STI})}]_{q\to
0}\,,$ (14)

is the theoretical predictions for the same form factor when the SM is turned
off.

Most importantly, as has been established recently in Aguilar  _et al._
(2023a), the difference between $\mathit{L}_{{sg}}(r^{2})$ and
$\mathit{L}_{{sg}}^{\\!\star}(r^{2})$ corresponds to a clearly nonvanishing
${\mathbb{C}}(r^{2})$, whose shape is absolutely compatible with that obtained
in (i).

Figure 3:  Left panel: Quenched lattice data for $\mathit{L}_{{sg}}(r^{2})$
from Aguilar  _et al._ (2021a) (points + black continuous curve) compared to
the $\mathit{L}_{{sg}}^{\\!\star}(r^{2})$ (green dashed) determined in Aguilar
_et al._ (2023a). The green band represents the error in
$\mathit{L}_{{sg}}^{\\!\star}(r^{2})$, estimated as explained in detail in
Aguilar  _et al._ (2023a). Right panel: The displacement function,
${\mathbb{C}}(r^{2})$, obtained from the results on the left panel through Eq.
(13) (points + black continuous curve), compared to the BS amplitude of
Aguilar  _et al._ (2022) (purple line and corresponding error band).

As we will see in the next section, the SM produces a displacement analogous
to that of Eq. (13) in the form factors of the quark-gluon vertex.

${\mathbb{C}}(r^{2})$
---
Defining equation | $\displaystyle{\lim_{q\to 0}}\frac{V_{1}(q,r,k)}{q^{2}}=\frac{2(q\cdot r)}{q^{2}}{\mathbb{C}}(r^{2})$ | ${\mathbb{C}}(r^{2})=\int\\!\small{d}^{{\scriptscriptstyle 4}}\\!\small{\ell}\,{\cal K}(r,\ell){\mathbb{C}}(\ell^{2})$ | $\mathit{L}_{{sg}}(r^{2})=\mathit{L}_{{sg}}^{\\!\star}(r^{2})+{\mathbb{C}}(r^{2})$
Terminology |  Momentum-dependent residue or residue function | BS amplitude | Displacement function
Table 1:  Equivalent definitions of the function ${\mathbb{C}}(r^{2})$ and the
corresponding terminology.

We close this introductory section by pointing out that, due to its multiple
roles, the function ${\mathbb{C}}(r^{2})$ receives a variety of names, which,
depending on the situation, emphasize different physical aspects; the
terminology employed and its origin is summarized in Table 1. Completely
analogous terminology is employed for the corresponding quantities associated
with the quark-gluon vertex, introduced in the next section.

##  III Quark-gluon vertex and displaced soft-gluon limit

The previous considerations have been carried out in the realm of pure Yang-
Mills theories, or quarkless QCD. We now turn to the main part of this study,
and consider a real-world version of QCD, with two light quark flavors.

In this case, the main modification with respect to the Yang-Mills case is the
inclusion of the fully-dressed quark-gluon vertex,
${\rm{I}}\\!\Gamma^{a}_{\alpha}(q,p_{2},-p_{1})$, which may be cast in the
form

${\rm{I}}\\!\Gamma_{\alpha}^{a}(q,p_{2},-p_{1})=ig\,\frac{\lambda^{a}}{2}\,{\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})\,,$
(15)

where the $\lambda^{a}$ denote the standard Gell-Mann matrices, and the four-
momenta satisfy the relation $q+p_{2}-p_{1}=0$.

At tree-level, ${\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})$ acquires the
standard expression

$\Gamma_{\\!0}^{\alpha}(q,p_{2},-p_{1})=\gamma^{\alpha}\,.$ (16)

The activation of the SM generates longitudinal poles in the only channel that
carries a Lorentz index, namely the $q$-channel. Then, as shown in the bottom
row of Fig. 2, the quark-gluon vertex is composed by two distinct pieces,

${\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})=\Gamma_{\alpha}(q,p_{2},-p_{1})+V_{\alpha}(q,p_{2},-p_{1})\,,$
(17)

where $\Gamma_{\alpha}(q,p_{2},-p_{1})$ is the pole-free part, and

$V_{\alpha}(q,p_{2},-p_{1})=\frac{q_{\alpha}}{q^{2}}Q(q,p_{2},-p_{1})\,,$ (18)

is the pole term.

The amplitude $Q(q,p_{2},-p_{1})$ may be decomposed in a standard Dirac basis
according to

$Q(q,p_{2},-p_{1})=Q_{1}\,\mathbb{I}+Q_{2}\,\not{p}_{2}+Q_{3}\,\not{p}_{1}+Q_{4}\,\tilde{\sigma}_{\mu\nu}\,p_{2}^{\mu}\,p_{1}^{\nu}\,,$
(19)

where $\mathbb{I}$ is the $4\times 4$ identity matrix in the Dirac space,
$\tilde{\sigma}^{\mu\nu}=1/2[\gamma^{\mu},\gamma^{\nu}]$, and we use the
short-hand notation $Q_{i}\equiv Q_{i}(q^{2},p_{2}^{2},p_{1}^{2})$.

###  III.1 Schwinger mechanism turned off

Let us suppose that the SM is not active, such that
$V_{\alpha}(q,p_{2},-p_{1})=0$; then, from Eq. (17) we have that
${\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})=\Gamma_{\alpha}(q,p_{2},-p_{1})$,
i.e., the full vertex is pole-free. Then, consider the STI triggered when the
quark-gluon vertex $\Gamma_{\alpha}(q,p_{2},-p_{1})$ is contracted by the
gluon momentum $q^{\alpha}$, namely Roberts and Williams (1994); Cloet and
Roberts (2014)

$q^{\alpha}\Gamma_{\alpha}(q,p_{2},-p_{1})=F(q^{2})[S^{-1}(p_{1})H(q,p_{2},-p_{1})-\overline{H}(-q,p_{1},-p_{2})\,S^{-1}(p_{2})]\,,$
(20)

where $S^{ab}(p)=i\delta^{ab}S(p)$ is the quark propagator Natale and
Rodrigues da Silva (1997); Fleischer  _et al._ (1999); Franco and Lubicz
(1998); Maris and Roberts (2003); Bhagwat  _et al._ (2003); Bowman  _et al._
(2005); Furui and Nakajima (2006); Sauli  _et al._ (2007); Chang  _et al._
(2007); Bashir  _et al._ (2012); Oliveira  _et al._ (2019b); Chen  _et al._
(2021); Gao  _et al._ (2021)

$S^{-1}(p)=A(p^{2})\,\not{p}-B(p^{2})\,.$ (21)

In addition, $H(q,p_{2},-p_{1})$ is the quark-ghost scattering kernel [see
upper panel of Fig. 4], and $\overline{H}(-q,p_{1},-p_{2})$ its “conjugate”,
whose relation to $H$ is explained in detail in Aguilar  _et al._ (2017) [see
discussion around Eq. (2.5)]. Note that the color structure has been factored
out, setting $H^{a}=-g\,(\lambda^{a}/2)H$.

Next, we expand both sides of Eq. (20) around $q=0$, with $p_{1}=p_{2}=p$. It
is clear that the zeroth order term vanishes on both sides; then, the matching
of the coefficients multiplying the terms linear in $q$ yields the
corresponding WI.

We start by noting that the term linear in $q$ on the l.h.s of Eq. (20) is
simply $q^{\alpha}\Gamma_{\alpha}(0,p,-p)$, where $\Gamma_{\alpha}(0,p,-p)$ is
the soft-gluon limit of the quark-gluon vertex. The most general Lorentz
decomposition of the vertex $\Gamma_{\alpha}(0,p,-p)$ is given by

$\Gamma_{\alpha}(0,p,-p)=\Gamma_{1}(p^{2})\,\gamma_{\alpha}+4\Gamma_{2}(p^{2})\,\not{p}\,p_{\alpha}+2\Gamma_{3}(p^{2})\,p_{\alpha}+2\Gamma_{4}(p^{2})\,\tilde{\sigma}_{\alpha\nu}p^{\nu}\,,$
(22)

where $\Gamma_{i}(p^{2})$ are the form factors in the soft-gluon limit.
However, due to charge conjugation symmetry Kızılersü  _et al._ (2021),
$\Gamma_{4}(p^{2})$ vanishes identically, thus reducing the number of relevant
form factors down to three.

Figure 4: Top panel: Diagrammatic representation of the quark-ghost kernel,
$H^{a}(q,p_{2},-p_{1})$. Bottom panel: The one-loop dressed approximation of
$H^{a}(q,p_{2},-p_{1})$, employed for the determination of the $K_{i}(p^{2})$
in Sec. V.

In the Landau gauge, the expansion of the r.h.s. of Eq. (20) around $q=0$ is
facilitated by setting Aguilar  _et al._ (2014)

$\displaystyle H(q,p_{2},-p_{1})=$
$\displaystyle\,Z_{H}+q^{\alpha}K_{\alpha}(q,p_{2},-p_{1})\,,$
$\displaystyle\overline{H}(-q,p_{1},-p_{2})=$
$\displaystyle\,Z_{H}-q^{\alpha}\overline{K}_{\alpha}(-q,p_{1},-p_{2})\,,$
(23)

where $Z_{H}$ is the quark-ghost kernel renormalization constant, which is
finite in Landau gauge Taylor (1971). The forms of the $K_{\alpha}$ and
$\overline{K}_{\alpha}$ have been given in Aguilar  _et al._ (2014). Thus, the
first two terms in the Taylor expansion of these kernels are given by

$\displaystyle H(0,p,-p)=$
$\displaystyle\,Z_{H}+q^{\alpha}K_{\alpha}(0,p,-p)\,,$
$\displaystyle\overline{H}(0,p,-p)=$
$\displaystyle\,Z_{H}-q^{\alpha}\overline{K}_{\alpha}(0,p,-p)\,,$ (24)

where, employing the same basis as in Eq. (22), we find

$\displaystyle K_{\mu}(0,p,-p)=$
$\displaystyle\,K_{1}(p^{2})\,\gamma_{\mu}+4K_{2}(p^{2})\,\not{p}\,p_{\mu}+2K_{3}(p^{2})\,p_{\mu}+2K_{4}(p^{2})\,\tilde{\sigma}_{\mu\nu}p^{\nu}\,,$
$\displaystyle\overline{K}_{\mu}(0,p,-p)=$
$\displaystyle\,K_{1}(p^{2})\,\gamma_{\mu}+4K_{2}(p^{2})\,\not{p}\,p_{\mu}+2K_{3}(p^{2})\,p_{\mu}-2K_{4}(p^{2})\,\tilde{\sigma}_{\mu\nu}p^{\nu}\,,$
(25)

with $K_{i}(p^{2})$ being the form factors. In addition, we use the following
expansion for the inverse quark propagator,

$S^{-1}(p_{1})\\!=\\!S^{-1}(p)\,+\,q^{\alpha}\bigg{[}A(p^{2})\gamma_{\alpha}+2A^{\prime}(p^{2})\,\not{p}\,p_{\alpha}-2B^{\prime}(p^{2})p_{\alpha}\bigg{]}\,,$
(26)

where the short-hand notation $f^{\prime}(p^{2}):={df(p^{2})}/{dp^{2}}$ is
used.

Then, the matching of the terms on the two sides of Eq. (20) expresses the
form factors $\Gamma_{i}(p^{2})$ in terms of the functions
$\\{A,B,F,K_{i}\\}$. Converting the results to Euclidean space by means of
standard transformation rules, and setting
$\lambda_{i}^{\star}(p^{2}):=\Gamma_{i}^{E}(p^{2}_{{\scriptscriptstyle E}})$,
we obtain

$\displaystyle\lambda_{1}^{\star}(p^{2})=$
$\displaystyle\,F(0)\left\\{\left[Z_{H}+4p^{2}K_{4}(p^{2})\right]A(p^{2})-2K_{1}(p^{2})B(p^{2})\right\\}\,,$
$\displaystyle\lambda_{2}^{\star}(p^{2})=$
$\displaystyle\,F(0)\left\\{-\frac{1}{2}A^{\prime}(p^{2})Z_{H}+\left[K_{3}(p^{2})+K_{4}(p^{2})\right]A(p^{2})+2K_{2}(p^{2})B(p^{2})\right\\}\,,$
$\displaystyle\lambda_{3}^{\star}(p^{2})=$
$\displaystyle\,F(0)\left\\{Z_{H}B^{\prime}(p^{2})+\left[K_{1}(p^{2})+4p^{2}K_{2}(p^{2})\right]A(p^{2})-2K_{3}(p^{2})B(p^{2})\right\\}\,.$
(27)

Note that Eq. (27) coincides with the Euclidean form of Eq. (4.3) in Aguilar
_et al._ (2014) up to redefinitions444The $K_{i}$ of Aguilar  _et al._ (2014)
can be obtained from the present ones by the redefinitions $K_{1}\to K_{1}$,
$K_{2}\to K_{2}/4$, $K_{3}\to K_{3}/2$ and $K_{4}\to K_{2}/2$.. Also, note
that Aguilar  _et al._ (2014) employed explicitly the Taylor renormalization
scheme, where $Z_{H}=1$ Taylor (1971).

###  III.2 Schwinger mechanism turned on

Let us now restore the full content of Eq. (17), i.e., allow for the presence
of massless poles in the fundamental vertices of the theory. With the
Schwinger mechanism activated, the form of the STI remains intact Smit (1974);
Aguilar  _et al._ (2016a): the divergence of the full vertex, now
${\rm{I}}\\!\Gamma_{\alpha}=\Gamma_{\alpha}+V_{\alpha}$, is still given by the
r.h.s. of Eq. (20). Since $q^{\alpha}V_{\alpha}=Q$, the STI satisfied by the
pole-free part of the vertex (the one simulated on the lattice) is displaced
with respect to Eq. (20), namely

$q^{\alpha}\Gamma_{\alpha}(q,p_{2},-p_{1})=F(q^{2})[S^{-1}(p_{1})H(q,p_{2},-p_{1})-\overline{H}(-q,p_{1},-p_{2})\,S^{-1}(p_{2})]-\underbrace{Q(q,p_{2},-p_{1})}_{\rm
displacement}.$ (28)

Consequently, the form factors of the vertex $\Gamma_{\alpha}$ in the presence
of the Schwinger mechanism, to be denoted by $\lambda_{i}$, will be also
displaced with respect to the $\lambda_{i}^{\star}$ given in Eq. (27). To
determine this displacement in detail, the construction leading to Eq. (27)
must be supplemented by the linear contributions arising from the Taylor
expansion of $Q(q,p_{2},-p_{1})$ around $q=0$.

To that end, note that, since at $q=0$ all other terms in Eq. (28) vanish, the
condition $Q(0,p,-p)=0$ must be fulfilled. Substituting $q=0$ into Eq. (19),
and setting $Q_{i}(0,p^{2},p^{2}):=Q_{i}(p^{2})$, we obtain the constraints

$Q_{1}(p^{2})=Q_{2}(p^{2})+Q_{3}(p^{2})=0\,.$ (29)

This argument leaves $Q_{4}(p^{2})$ undetermined; however, charge conjugation
symmetry forces $Q_{4}(p^{2})$ to vanish, $Q_{4}(p^{2})=0$.

Thus, the Taylor expansion of $Q(q,p_{2},-p_{1})$ is given by

$\lim_{q\to 0}Q(q,p_{2},-p_{1})=q^{\alpha}\left[\frac{\partial
Q(q,p_{2},-p_{1})}{\partial q^{\alpha}}\right]_{q=0}+{\cal O}(q^{2})\,,$ (30)

and from Eq. (19) we find

$\left[\frac{\partial Q(q,p_{2},-p_{1})}{\partial
q^{\alpha}}\right]_{q=0}=2{\mathbb{Q}}_{1}(p^{2})\,p_{\alpha}+2{\mathbb{Q}}_{2+3}(p^{2})\,\not{p}\,p_{\alpha}+Q_{3}(p^{2})\,\gamma_{\alpha}\,,$
(31)

where we have introduced

$\displaystyle{\mathbb{Q}}_{i}(p^{2}):=$ $\displaystyle[\partial
Q(q,p_{2},-p_{1})/\partial p_{1}^{2}]_{q=0}\,,\,\qquad i=1,2,3,$
$\displaystyle{\mathbb{Q}}_{2+3}(p^{2}):=$
$\displaystyle{\mathbb{Q}}_{2}(p^{2})+{\mathbb{Q}}_{3}(p^{2})\,.$ (32)

Then, combining Eqs. (27) and (30), we may easily determine from Eq. (28) the
displacement of the vertex form factors when the SM is active, namely

$\displaystyle\lambda_{1}(p^{2})=$
$\displaystyle\,\lambda_{1}^{\star}(p^{2})\,-\,Q_{3}(p^{2})\,,$
$\displaystyle\lambda_{2}(p^{2})=$
$\displaystyle\,\lambda_{2}^{\star}(p^{2})\,+\,\frac{1}{2}{\mathbb{Q}}_{2+3}(p^{2})\,,$
$\displaystyle\lambda_{3}(p^{2})=$
$\displaystyle\,\lambda_{3}^{\star}(p^{2})\,+\,{\mathbb{Q}}_{1}(p^{2})\,.$
(33)

The above relations are the direct analogue of Eqs. (13) and (14) for the case
of the quark-gluon vertex; $Q_{3}(p^{2})$, ${\mathbb{Q}}_{2+3}(p^{2})$, and
${\mathbb{Q}}_{1}(p^{2})$ are the corresponding displacement functions.

We emphasize that the $\lambda_{i}(p^{2})$ capture the full dynamical content
of the theory, in the sense that they take into account the action of the SM.
Therefore, if the SM is realized as described here, it is the
$\lambda_{i}(p^{2})$, and not the $\lambda_{i}^{\star}(p^{2})$, that should
coincide with the results obtained for these form factors from lattice QCD.
Let us now suppose that the $\lambda_{i}^{\star}(p^{2})$ were determined
through Eq. (27), using lattice ingredients for all (or most) intermediate
calculations. Then, if a statistically significant discrepancy is found
between the $\lambda_{i}^{\star}(p^{2})$ and the lattice $\lambda_{i}(p^{2})$,
it is natural to attribute it to the presence of the displacement functions.
Such an observation becomes even more compelling if the observed signal
displays a momentum dependence compatible with the BS prediction for the
corresponding amplitudes/displacement functions, as happens in the case of
pure Yang-Mills Aguilar  _et al._ (2023a) (see Fig. 3). As we demonstrate in
the rest of this article, this is indeed what happens also in the case of QCD
with $N_{f}=2$.

##  IV Dynamical determination of the displacement functions

In order to determine the numerical importance of the displacements exhibited
in Eq. (33), in this section we compute the displacement functions from the BS
equations that they obey.

###  IV.1 System of Bethe-Salpeter equations

The starting point of this analysis is the coupled system of SDEs satisfied by
the vertices ${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)$ and
${\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})$, shown in Fig. 5. Note that we
opt for the version of the SDE obtained within the 3PI formalism; as a result,
the vertices carrying the momentum $q$ are fully dressed. The corresponding
four-particle kernels ${\cal K}_{ij}$ are appropriately modified to avoid
overcounting; as may be deduced from Fig. 6, for the purposes of this
computation the “one-particle exchange” approximations of these kernels will
be employed.

Figure 5:  The SDEs for the three-gluon (top row) and quark-gluon (bottom row)
vertices. The white (colored) circles denote fully-dressed propagators
(vertices), while the gray ellipses denote four-particle kernels.

To obtain the relevant dynamical equations, we substitute for the vertices
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)$ and
${\rm{I}}\\!\Gamma_{\alpha}(q,p_{2},-p_{1})$ appearing on either side of the
SDE system (red and orange circles) the expressions given in Eqs. (4) and
(17), respectively. This introduces the pole parts on both sides of the
system, and, after contraction of the first SDE (top row) by ${\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,k)$, one may use the expressions
given by Eqs. (5) and (18). Then, as the limit $q\to 0$ is taken on both sides
of the SDEs, the pole terms dominate over the regular parts, and, after
triggering Eqs. (6) and (30), the displacement functions emerge as the leading
contributions. Finally, the proper identification of the various tensorial
structures on both sides, and subsequent matching of their co-factors, gives
rise to a BS system that may be schematically represented as

${\cal A}_{m}=\sum_{n=1}^{4}\int\\!\\!d^{\,4}\ell\,{\cal M}_{mn}\,{\cal
A}_{n}\quad;\quad m=1,2,3,4\,,$ (34)

where the four-entry quantity ${\cal
A}=\big{[}{\mathbb{C}},Q_{3},{\mathbb{Q}}_{2+3},{\mathbb{Q}}_{1}\big{]}$ was
introduced, and the functions ${\cal M}_{mn}$ contain all remaining
ingredients; their explicit form depends on the approximations employed for
the kernels ${\cal K}_{ij}$.

Figure 6:  The system of coupled BS equations that determines the amplitudes
${\mathbb{C}}(k^{2})$ and $Q_{3}(p^{2})$, in the one-particle exchange
approximation. In the first line, a summation of the quark flavors is implied,
while the factor of 2 accounts for the two orientations of the quark loop.

In what follows we approximate the kernels ${\cal K}_{ij}$ by their one-
particle exchange form, shown in Fig. 6. In addition, we focus exclusively on
the form factor $\lambda_{1}$, associated with the tree-level tensor
$\gamma_{\alpha}$, and determine the displacement function $Q_{3}$ appearing
in the first relation of Eq. (33). Therefore, we restrict our analysis to the
reduced system comprised by ${\mathbb{C}}(r^{2})$ and $Q_{3}(p^{2})$ only, as
shown in Fig. 6. We have confirmed that the omission of the functions
${\mathbb{Q}}_{2+3}(p^{2})$ and ${\mathbb{Q}}_{1}(p^{2})$ has practically no
impact on ${\mathbb{C}}(r^{2})$ and $Q_{3}(p^{2})$.

In order to proceed, we need to furnish inputs for the Landau-gauge
propagators and transversely projected vertices appearing in the system of
Fig. 6. For the gluon, ghost, and quark propagators, we use appropriate fits
to results obtained from lattice QCD, as described in App. A. Furthermore, we
assume that the transversely projected three-gluon and quark-gluon vertices,

$\displaystyle\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{\alpha\mu\nu}(q,r,k):=$
$\displaystyle\,P_{\alpha}^{\alpha^{\prime}}(q)P_{\mu}^{\mu^{\prime}}(r)P_{\nu}^{\nu^{\prime}}(k){\rm{I}}\\!\Gamma_{\alpha^{\prime}\mu^{\prime}\nu^{\prime}}(q,r,k)\,,$
$\displaystyle\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{\alpha}(q,p_{2},-p_{1}):=$
$\displaystyle\,P^{\alpha^{\prime}}_{\alpha}(q){\rm{I}}\\!\Gamma_{\alpha^{\prime}}(q,p_{2},-p_{1})\,,$
(35)

respectively, can be approximated by

$\displaystyle\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}^{\alpha\mu\nu}(q,r,k)=$
$\displaystyle\,\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{0}^{\alpha\mu\nu}(q,r,k)\mathit{L}_{{sg}}(s^{2})\,,\qquad
s^{2}=\frac{1}{2}(q^{2}+r^{2}+k^{2})\,,$
$\displaystyle\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}^{\alpha}(q,p_{2},-p_{1})=$
$\displaystyle\,\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{0}^{\alpha}(q,p_{2},-p_{1})\lambda_{1}(\bar{s}^{2})\,,\qquad\bar{s}^{2}=\frac{1}{2}(q^{2}+p_{1}^{2}+p_{2}^{2})\,,$
(36)

where
$\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{0}^{\alpha\mu\nu}$
and $\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{0}^{\alpha}$
are the tree-level forms of
$\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{\alpha\mu\nu}$ and
$\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{\alpha}$, and are
obtained by transversely projecting the $\Gamma_{\\!0}^{\alpha\mu\nu}$ and
$\Gamma_{\\!0}^{\alpha}$ of Eqs. (12) and (16), respectively. For the three-
gluon vertex, Eq. (36) has been validated by numerous studies in pure Yang-
Mills theory Eichmann  _et al._ (2014); Blum  _et al._ (2014); Huber (2020a);
Pinto-Gómez  _et al._ (2023); Aguilar  _et al._ (2023c), and we assume it to
hold as well for $N_{f}=2$. In the case of the quark-gluon vertex, the
particular form in Eq. (36) has not been tested explicitly in the literature.
However, previous works have already shown that the classical form factor of
$\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}_{\alpha}(q,p_{2},-p_{1})$
is the largest in magnitude Skullerud and Kizilersu (2002); Skullerud  _et
al._ (2003, 2005); Lin (2006); Alkofer  _et al._ (2009); Aguilar  _et al._
(2014); Mitter  _et al._ (2015); Cyrol  _et al._ (2018); Gao  _et al._ (2021);
Kızılersü  _et al._ (2021), and that its angular dependence is rather weak
Kizilersu  _et al._ (2007); Blum (2017); Blum  _et al._ (2015); Cyrol  _et
al._ (2018); together, these observations motivate the use of Eq. (36) for
$\overline{{\rm{I}}\\!\Gamma}\vphantom{{\rm{I}}\\!\Gamma}^{\alpha}(q,p_{2},-p_{1})$.

In order to unify the description of the results, we rename $r^{2}\to p^{2}$.
Next, we pass to Euclidean momenta, i.e., $p^{2}\to-p^{2}_{{\scriptscriptstyle
E}}$, with $p^{2}_{{\scriptscriptstyle E}}>0$, introduce the integral measure
in hyperspherical coordinates

$\int\\!d^{4}\\!\ell_{{\scriptscriptstyle
E}}=\frac{1}{(2\pi)^{3}}\int_{y,\theta}\quad;\qquad\int_{y,\theta}:=\int_{0}^{\infty}\\!\\!dy\,y\int_{0}^{\pi}\\!\\!d\theta\,s^{2}_{\theta}\,,$
(37)

and employ standard transformation properties for all Green’s functions and
form factors involved (see, e.g., Eqs. (5.13) to (5.16) of Aguilar  _et al._
(2014)); note, in particular, that ${\mathbb{C}}^{{\scriptscriptstyle
E}}(p^{2}_{{\scriptscriptstyle E}})=-{\mathbb{C}}(-p^{2}_{{\scriptscriptstyle
E}})$ Aguilar  _et al._ (2022) and $Q^{{\scriptscriptstyle
E}}_{3}(p^{2}_{{\scriptscriptstyle E}})=Q_{3}(-p^{2}_{{\scriptscriptstyle
E}})$.

Defining

$\displaystyle x=p^{2}\,,\quad y=\ell^{2}\,,\quad
p\cdot\ell=\sqrt{xy}c_{\theta}\,,$ $\displaystyle
z=(p+\ell)^{2}=x+y+2\sqrt{xy}c_{\theta}\,,$ $\displaystyle
u=(p-\ell)^{2}=x+y-2\sqrt{xy}c_{\theta}\,,$ (38)

where the index “E” has been suppressed throughout, we arrive at the following
system of coupled BS equations

$\displaystyle{\mathbb{C}}(x)=$
$\displaystyle\,-\frac{\alpha_{s}}{\pi^{2}}\int_{y,\theta}\\!\\!f_{11}\,{\mathbb{C}}(y)\,+\,\frac{\alpha_{s}}{6\pi^{2}}\int_{y,\theta}\\!\\!f_{12}\,Q_{3}(y)\,,$
$\displaystyle Q_{3}(x)=$
$\displaystyle\,\frac{3\alpha_{s}}{4\pi^{2}}\int_{y,\theta}\\!\\!f_{21}\,{\mathbb{C}}(y)\,+\,\frac{\alpha_{s}}{12\pi^{2}}\int_{y,\theta}\\!\\!f_{22}\,Q_{3}(y)\,.$
(39)

The functions $f_{ij}(x,y,\theta)$ are given by

$\displaystyle f_{11}=$
$\displaystyle\,\frac{s_{\theta}^{2}c_{\theta}}{z}\sqrt{\frac{y}{x}}\left[3(x^{2}+y^{2})+6c_{\theta}\sqrt{xy}(x+y)+xy(c_{\theta}^{2}+8)\right]\Delta(z)\Delta^{2}(y)\mathit{L}_{{sg}}^{2}(s^{2})\,,$
$\displaystyle f_{12}=$
$\displaystyle\,4\sqrt{\frac{y}{x}}c_{\theta}\left\\{3A_{y}B_{y}B_{z}+\big{[}(3-2s^{2}_{\theta})y+3\sqrt{xy}c_{\theta}\big{]}A^{2}_{y}A_{z}\right\\}\frac{\lambda_{1}^{2}(s^{2})}{(yA^{2}_{y}+B^{2}_{y})^{2}(zA^{2}_{z}+B^{2}_{z})}$
$\displaystyle\,\hskip
14.22636pt-6\bigg{(}1+\sqrt{\frac{y}{x}}c_{\theta}\bigg{)}\frac{A_{z}\lambda_{1}^{2}(s^{2})}{(yA^{2}_{y}+B^{2}_{y})(zA^{2}_{z}+B^{2}_{z})}\,,$
$\displaystyle f_{21}=$
$\displaystyle\,-2ys^{2}_{\theta}\bigg{(}1+\frac{2}{3}\sqrt{\frac{x}{y}}c_{\theta}\bigg{)}\frac{A_{z}}{zA^{2}_{z}+B^{2}_{z}}\,\Delta^{2}(y)\lambda_{1}^{2}(s^{2})\,,$
$\displaystyle f_{22}=$
$\displaystyle\,\left\\{\frac{2}{3}y\bigg{[}1+\frac{2}{u}(y-\sqrt{xy}c_{\theta})\bigg{]}\frac{s^{2}_{\theta}A^{2}_{y}}{(yA^{2}_{y}+B^{2}_{y})}-1-\frac{2}{3}\frac{y}{u}s^{2}_{\theta}\right\\}\frac{\Delta(u)\lambda_{1}^{2}({\bar{s}}^{2})}{(yA^{2}_{y}+B^{2}_{y})}\,,$
(40)

where $s^{2}=x+y+\sqrt{xy}c_{\theta}$,
${\bar{s}}^{2}=x+y-\sqrt{xy}c_{\theta}$, and the subscripts $y$ or $z$ in the
functions $A$ and $B$ indicate their momentum dependence i.e., $A_{z}:=A(z)$,
etc.

###  IV.2 Inputs

In order to proceed with the solution of the system given by Eq. (39), we need
to specify the form of the various components entering in the functions
$f_{ij}(x,y,\theta)$, and in particular $A$, $B$, $\Delta$,
$\mathit{L}_{{sg}}$, and $\lambda_{1}$. In addition, the ghost dressing
function $F$ is needed for implementing the required conversions between
renormalization schemes, as explained below.

Note in particular the following important points:

(i) Throughout this work, we will use as external inputs a series of fits
given in App. A to the $N_{f}=2$ lattice data for the gluon and ghost
propagators from Ayala  _et al._ (2012); Binosi  _et al._ (2017b), and for
$A(p^{2})$, $B(p^{2})$ and $\lambda_{1}(p^{2})$ from Oliveira  _et al._
(2019b); Kızılersü  _et al._ (2021). Specifically, for the quark functions, we
employ the setups denominated “L08” and “L07” in Table I of Kızılersü  _et
al._ (2021), which have large statistics and small current quark mass,
$m_{q}$, and pion mass, $m_{\pi}$, namely $(m_{q},m_{\pi})=(6.2,280)$ MeV and
$(8,295)$ MeV, respectively.

(ii) The lattice data of Oliveira  _et al._ (2019b); Kızılersü  _et al._
(2021) for the functions $A(p^{2})$ and $\lambda_{1}(p^{2})$ display visible
artifacts in the ultraviolet, where the known perturbative behavior of these
functions is not accurately reproduced. This issue, which has been discussed
in various studies Skullerud  _et al._ (2001, 2005, 2003); Lin (2006);
Sternbeck  _et al._ (2017); Kizilersu  _et al._ (2007); Oliveira  _et al._
(2019b); Kızılersü  _et al._ (2021), is ameliorated by the use of the so-
called “tree-level correction” Skullerud  _et al._ (2001); Oliveira  _et al._
(2019b); Kızılersü  _et al._ (2021), but even then is not completely cured. In
fact, signs of these artifacts are already present at $p=3$ GeV, where the L07
$\lambda_{1}(p^{2})$ starts to increase [see bottom left panel of Fig. 10],
rather than continuing decreasing as predicted by perturbation theory Pascual
and Tarrach (1984); Muta (1987); Davydychev  _et al._ (2001). For this reason,
we discard the data for these functions for momenta above $2.5$ GeV. Since we
employ fitting functions that reproduce the one-loop resumed perturbative
result for each of the functions at large momenta (see App. A), the
ultraviolet behavior of our inputs is under control.

(iii) We next turn to the $\mathit{L}_{{sg}}(p^{2})$. Unfortunately, the only
available unquenched lattice results for $\mathit{L}_{{sg}}(r^{2})$ are not
for $N_{f}=2$ but rather for $N_{f}=2+1$ (two light quarks with current mass
$1.3$ MeV, and a heavier one, with current mass $63$ MeV) Aguilar  _et al._
(2020). However, at least in the case of the gluon propagator Ayala  _et al._
(2012); Aguilar  _et al._ (2020), the difference between $N_{f}=2$ to
$N_{f}=2+1$ is rather mild. It is therefore reasonable to use for the
$\mathit{L}_{{sg}}(p^{2})$ in the present analysis the $N_{f}=2+1$ data from
Aguilar  _et al._ (2020). We hasten to emphasize that this last approximation
only affects the pole amplitudes determined through the BS equations, whereas
the determination of $Q_{3}(p^{2})$, carried out in the next section through
the analysis of the WI displacement, does not depend on
$\mathit{L}_{{sg}}(p^{2})$.

(iv) In general, the aforementioned lattice data are renormalized in different
schemes; therefore, their self-consistent use hinges on their careful
conversion into a common renormalization scheme. To that end, we adopt a
particular variation of momentum subtraction (MOM) scheme, the so-called
$\widetilde{\text{MOM}}$ scheme Skullerud and Kizilersu (2002), defined by the
prescription

$\Delta^{-1}(\mu^{2})=\mu^{2}\,,\quad F(\mu^{2})=1\,,\quad
A(\mu^{2})=1\,,\quad\lambda_{1}(\mu^{2})=1\,,$ (41)

where we take $\mu=2$ GeV. The procedure employed for consistently adjusting
the results between renormalization schemes is described in detail in App. B.

Figure 7:  The green continuous curve on the left panel and the black
continuous curve (with blue band) on the right are the solutions obtained from
Eq. (39) using unquenched ($N_{f}=2$) lattice results as inputs. The input
used for $\lambda_{1}(s^{2})$ is taken from the L08 setup of Kızılersü  _et
al._ (2021); the results using LO7 inputs are nearly identical. The purple
dot-dashed curve on the left panel is the corresponding solution when quenched
lattice inputs are employed (pure Yang-Mills case), and is displayed for the
purpose of comparison.

###  IV.3 Results

With all external ingredients consistently renormalized in the
$\widetilde{\text{MOM}}$ scheme, we are in position to solve the BS equations
of Eq. (39) numerically, converting it to an eigenvalue problem. There are two
important points to emphasize. First, as happened in the quenched case Aguilar
_et al._ (2022), the coupling strength required for obtaining a solution
exceeds the $\widetilde{\text{MOM}}$ value; specifically, we have that
$\alpha^{\rm MOM}_{s}=0.47$ [see Eq. (66)], while the required value for
solving the system is $\alpha^{\rm BS}_{s}=1.17$ . The eigenvalue is known to
depend strongly on the details of the truncation used for the kernels ${\cal
K}_{ij}$, but affects only slightly the form of the solutions obtained.
Second, since the system of BS equations in Eq. (39) is linear and
homogeneous, its overall scale is undetermined: the multiplication of a given
solution by an arbitrary constant yields another solution.

In Fig. 7 we show the resulting ${\mathbb{C}}(p^{2})$ (left) and
$Q_{3}(p^{2})$ (right) as blue and black curves, respectively. This particular
solution has its scale set by fitting an arbitrary solution of the BS
equations to the result of the WI displacement, to be described in the next
section. The bands around the curves result from the standard error in
determining this scale.

As we can see in Fig. 7, for $p<1$ GeV the ${\mathbb{C}}(p^{2})$ and
$Q_{3}(p^{2})$ have similar magnitudes, (but opposite signs), suggesting that
both amplitudes might have similar weights in the BS dynamics. However, the
terms containing ${\mathbb{C}}(p^{2})$ in Eq. (39) have much larger
prefactors. Moreover, we can see from Fig. 7 that ${\mathbb{C}}(p^{2})$ has
support over a much longer interval of momenta than $Q_{3}(p^{2})$, which
decreases rapidly for $p>1$ GeV. As a result, ${\mathbb{C}}(p^{2})$ dominates
the dynamics of the BS system. In fact, the present solution for
${\mathbb{C}}(p^{2})$ is rather similar to quenched results obtained
previously Aguilar  _et al._ (2012a, 2018, 2022), shown as the purple dot-
dashed curve in the left panel of Fig. 7. Finally, we point out that if
${\mathbb{C}}(p^{2})$ is set to zero in Eq. (39), the resulting BS equation
for $Q_{3}(p^{2})$ has no solutions for $\alpha_{s}>0$.

Figure 8: Left panel: Lattice data for the gluon propagator for $N_{f}=0$
(brown diamonds) Aguilar  _et al._ (2021b) and $N_{f}=2$ (orange circles)
Ayala  _et al._ (2012); Binosi  _et al._ (2017b) . The black dot-dashed and
blue continuous curves represent the fits given in Eq. (C11) of Aguilar  _et
al._ (2022) and Eq. (50), respectively. In the inset we show the corresponding
gluon dressing functions. Right panel: The kernels ${\cal
K}_{\mathbb{C}}^{N_{f}=0}(k^{2})$ (blue dot-dashed), ${\cal
K}_{\mathbb{C}}^{N_{f}=2}(k^{2})$ (orange continuous), and ${\cal
K}_{Q_{3}}^{N_{f}=2}(k^{2})$ (purple dotted). The inset shows ${\cal
K}_{{\mathbb{C}}}^{N_{f}=2}(k^{2})+{\cal K}_{Q_{3}}^{N_{f}=2}(k^{2})$ as a red
continuous curve, compared to ${\cal K}_{\mathbb{C}}^{N_{f}=0}(k^{2})$. The
error bands are obtained through appropriate propagation of the corresponding
errors in the ${\mathbb{C}}(p^{2})$ and $Q_{3}(p^{2})$ of Fig. 7.

The solutions obtained from the BS system indicate that the mass generating
mechanism known from the pure Yang-Mills case is only mildly affected by the
inclusion of two active quark flavors. To appreciate this point in some
detail, we turn to the gluon mass equation that emerges from the treatment of
the two diagrams shown in Fig. 1; the additional gluonic corrections stemming
from the two-loop diagrams will be neglected for the purposes of this
qualitative analysis

The procedure outlined in Papavassiliou (2022); Ferreira and Papavassiliou
(2023) for obtaining the contribution to the gluon mass from diagram $d_{1}$
by appealing to Eq. (6) can be straightforwardly extended to the case of
diagram $d_{2}$, where Eq. (30) must be employed. Keeping only the dominant
term $Q_{3}(p^{2})$, one obtains (Euclidean space)

$m^{2}=Z_{3}\int_{0}^{\infty}\\!\\!dy\,{\cal
K}_{\mathbb{C}}^{N_{f}}(y)\,+\,Z_{2}\int_{0}^{\infty}\\!\\!dy\,{\cal
K}_{Q_{3}}^{N_{f}}(y)\,,$ (42)

where the renormalization constants $Z_{2}$ and $Z_{3}$ are defined in Eq.
(56), and the kernels are given by

${\cal
K}_{\mathbb{C}}^{N_{f}}(y)=-\frac{9\alpha_{s}}{8\pi}y^{2}\Delta^{2}(y){\mathbb{C}}(y)\,,\qquad{\cal
K}_{Q_{3}}^{N_{f}}(y)=\frac{\alpha_{s}N_{f}}{4\pi}\left[\frac{y\left(yA_{y}^{2}+2B^{2}_{y}\right)Q_{3}(y)}{\left(yA_{y}^{2}+B^{2}_{y}\right)^{2}}\right]\,.$
(43)

Note that the dependence of ${\cal K}_{\mathbb{C}}^{N_{f}}(y)$ on $N_{f}$ is
only implicit, due to the $N_{f}$ dependence of $\Delta(p^{2})$,
${\mathbb{C}}(p^{2})$ and $\alpha_{s}$. Instead, ${\cal K}_{Q_{3}}^{N_{f}}(y)$
depends explicitly on $N_{f}$, as well as implicitly.

As is known from lattice simulation Ayala  _et al._ (2012); Binosi  _et al._
(2017b), the saturation point of the gluon propagator with two active quarks
is lower with respect to the quenched case, as may be seen on the left panel
of Eq. (8); this is tantamount to saying that the gluon mass increases when
quarks are introduced. The kernels given in Eq. (43) seem to capture this
tendency: the inclusion of the quarks leads to kernels with additional support
throughout the relevant region of integration.

The origin of this relative enhancement is twofold. First, as shown on the
right panel of Fig. 8, ${\cal K}_{\mathbb{C}}^{N_{f}=2}(y)$ (orange
continuous) increases slightly with respect to ${\cal
K}_{\mathbb{C}}^{N_{f}=0}(y)$ (blue dot-dashed) (11% at the corresponding
peaks). This increase is the net outcome of two opposing tendencies: The
unquenched gluon propagators suppress ${\cal K}_{\mathbb{C}}^{N_{f}=2}(y)$
with respect to ${\cal K}_{\mathbb{C}}^{N_{f}=0}(y)$, while the unquenched
coupling enhances it Ayala  _et al._ (2012). Second, the kernel ${\cal
K}_{Q_{3}}^{N_{f}}(y)$ (purple dotted), which is absent when $N_{f}=0$,
provides an additional small positive contribution, as shown on the right
panel of Fig. 8. The overall effect may be roughly appreciated if we set
$Z_{3}=Z_{2}=1$ in Eq. (42) and then add the two kernels (red continuous); the
end result, compared to the quenched case (blue dot-dashed), is shown as an
inset on the right panel of Fig. 8.

This additional support of the unquenched kernel represents a relatively small
increment with respect to the quenched kernel; therefore it is reasonable to
expect a moderate increase in the gluon mass, compatible with the lattice
findings of Ayala  _et al._ (2012); Binosi  _et al._ (2017b). The detailed
calculation of the effect requires proper renormalization and the inclusion in
Eq. (42) of two-loop dressed loops, not shown here; however, this task lies
beyond the scope of the present work.

##  V Displacement function from lattice data

In this section we will test the first relation in Eq. (33), which expresses
the displacement associated with the classical form factor as

$Q_{3}(p^{2})=\lambda_{1}^{\star}(p^{2})-\lambda_{1}(p^{2})\,.$ (44)

To that end, we use the first relation in Eq. (27) to compute the WI
prediction, $\lambda_{1}^{\star}(p^{2})$, and then subtract from it the
lattice data of Kızılersü  _et al._ (2021) for $\lambda_{1}(p^{2})$.

The evaluation of the WI prediction is rather subtle, mainly due to the
presence of the factor $Z_{H}$, the renormalization constant of the quark-
ghost kernel. $Z_{H}$ is finite in the Landau gauge; in fact, there exists a
renormalization scheme, namely the Taylor scheme [see Eq. (59)], where
$Z_{H}=1$. However, this scheme is _not_ the same as the
$\widetilde{\text{MOM}}$ employed on the lattice determination of the Green’s
functions appearing in Eq. (27). Therefore, in App. B we determine the
appropriate values for the $\widetilde{\text{MOM}}$ scheme with $\mu=2$ GeV to
be $Z_{H}=1.120(8)$ and $Z_{H}=1.121(9)$, for the L08 and L07 lattice setups,
respectively.

In order to determine the form factors $K_{1}(p^{2})$ and $K_{4}(p^{2})$, we
evaluate the one-loop dressed diagram shown in Fig. 4, and use Eq. (23). The
full ghost-gluon vertex appearing in this diagram has the general form

${\rm{I}}\\!\Gamma_{\\!\\!c}^{\mu}(r,k,q)=B_{1}(r^{2},k^{2},q^{2})r^{\mu}+B_{2}(r^{2},k^{2},q^{2})q^{\mu}\,,$
(45)

where $r$, $k$, and $q$ stand for the momenta of the antighost, ghost, and
gluon, respectively; at tree-level, $B_{1}^{0}=1$ and $B_{2}^{0}=0$.

Then, it is relatively straightforward to obtain the results

$\displaystyle K_{1}(x)\\!=$
$\displaystyle\frac{\alpha_{s}Z_{H}}{4\pi^{2}}\int_{y,\theta}(2+c_{\theta}^{2})F(y)\Delta(y)B_{1}(y,0,y)\frac{B_{z}\lambda_{1}(s^{2})}{y\left(zA^{2}_{z}+B_{z}^{2}\right)}\,,$
$\displaystyle K_{4}(x)\\!=$
$\displaystyle-\frac{\alpha_{s}Z_{H}}{8\pi^{2}}\int_{y,\theta}\left(2+3\sqrt{\frac{y}{x}}c_{\theta}+c_{\theta}^{2}\right)\Delta(y)F(y)B_{1}(y,0,y)\frac{A_{z}\lambda_{1}(s^{2})}{y\left(zA^{2}_{z}+B_{z}^{2}\right)}\,,$
(46)

where $B_{1}(y,0,y)$ corresponds to the so-called “soft-ghost limit” of the
form factor $B_{1}$.

The numerical evaluation of the above expressions proceeds through the
substitution of $A$, $B$, $\Delta$, $\mathit{L}_{{sg}}$, and $\lambda_{1}$ by
fits to the lattice data, exactly as described in Sec. IV. For $\alpha_{s}$,
we use the $\widetilde{\text{MOM}}$ values determined in Eq. (66). Regarding
$B_{1}(y,0,y)$, since no lattice data are available for this quantity for
$N_{f}=2$, we determine it from the SDE that it satisfies [see App. C].

Figure 9:  The dimensionless combinations $p^{2}K_{4}(p^{2})A(p^{2})$ (blue
continuous) and $-K_{1}(p^{2})B(p^{2})$ (red dashed) determining
$\lambda_{1}^{\star}(p^{2})$, viz. first relation in Eq. (27). The left and
right panels correspond to the use of inputs from the lattice setups L08
(left) and L07 (right).

The dimensionless combinations $p^{2}K_{4}(p^{2})A(p^{2})$ and
$-K_{1}(p^{2})B(p^{2})$ appearing in Eq. (27) are shown in Fig. 9 as blue
continuous and red dashed lines, respectively. The two panels correspond to
the results obtained when the inputs used for $\lambda_{1}$ in Eq. (46)
originate from the lattice setup L08 (left) or L07 (right). The comparison of
the two panels reveals that $p^{2}K_{4}(p^{2})A(p^{2})$ is nearly identical
for both lattices, whereas the value of $-K_{1}(p^{2})B(p^{2})$ at the origin
is slightly larger in magnitude when the $\lambda_{1}$ from the L07 setup is
used.

Evidently, the errors associated with the inputs used to evaluate Eq. (46) get
propagated to the resulting $K_{i}(p^{2})$. It turns out that the largest
source of uncertainty is the form factor $\lambda_{1}$. In order to estimate
the error that $\lambda_{1}$ introduces to the $K_{i}(p^{2})$, we repeat the
numerical evaluation of Eq. (46) with

$\lambda_{1}(p^{2})\to\lambda_{1}(p^{2})\pm\delta\lambda_{1}(p^{2})\,,$ (47)

where $\delta\lambda_{1}(p^{2})$ is the $1\sigma$ spread in
$\lambda_{1}(p^{2})$, shown as green bands in the left panels of Fig. 10, for
each of the lattice setups. Then, the error obtained is combined in quadrature
with all other intrinsic errors stemming from the remaining ingredients,
mainly $A(p^{2})$ and $B(p^{2})$, to furnish our estimate for the total error
of the quantities $p^{2}K_{4}(p^{2})A(p^{2})$ and $-K_{1}(p^{2})B(p^{2})$.

Figure 10:  Top left: (i) the L08 lattice data (points) for
$\lambda_{1}(p^{2})$ Kızılersü  _et al._ (2021), the fit of Eq. (55) (purple
continuous), and the $1\sigma$ confidence band; (ii) the WI prediction
$\lambda_{1}^{\star}(p^{2})$ (green dot-dashed line and band). Top right: The
displacement function $Q_{3}(p^{2})$ (points) obtained as the difference
between the two curves shown on the left, as dictated by Eq. (44). The black
solid line and blue band indicate the prediction of the BS equation. The
bottom panels show the same quantities but for the L07 setup.

With the $K_{i}(p^{2})$ in hand, we use the first line of Eq. (27) to compute
$\lambda_{1}^{\star}(p^{2})$. The result is shown as a green dot-dashed line
for each of the lattice setups, L08 and L07, on the top left and bottom left
panels of Fig. 10, respectively. The green bands around these curves are our
error estimates, which combine in quadrature the errors in all of the
ingredients appearing in Eq. (27). In particular, we included in the total
error budget a $5\%$ error estimate in the value of $F(0)=2.39$ obtained by
extrapolating the lattice data of Ayala  _et al._ (2012); Binosi  _et al._
(2017b) to the origin, $p=0$.

As is clear from the left panels of Fig. 10, our results for
$\lambda_{1}^{\star}(p^{2})$ are considerably larger than the input
$\lambda_{1}(p^{2})$, indicating the existence of a nontrivial displacement,
$Q_{3}(p^{2})$. On the right panels of Fig. 10 we show as points the results
for $Q_{3}(p^{2})$ obtained from Eq. (33), for the L08 and L07 setups on the
top and bottom, respectively. The error bars in the points for $Q_{3}(p^{2})$
result from combining in quadrature the standard error of the lattice
$\lambda_{1}(p^{2})$ and the estimated error in $\lambda_{1}^{\star}(p^{2})$.

We next evaluate the statistical significance of the signal obtained for
$Q_{3}(p^{2})$, in comparison to the “null hypothesis”, namely the case when
$Q_{3}\to Q_{3}^{0}=0$, due to the the absence of the Schwinger mechanism. To
that end, we compute separately for each of the lattice setups the $\chi^{2}$
deviation from the null hypothesis as

$\chi^{2}=\sum_{i=1}^{n_{p}}\frac{[Q_{3}(p_{i}^{2})-Q_{3}^{0}(p_{i}^{2})]^{2}}{\epsilon_{Q_{3}(p_{i}^{2})}^{2}}\,,$
(48)

where the sum over the lattice points $i$, corresponding to momenta $p_{i}$,
extends until $p_{i}<2.5$ GeV; the reason for implementing this cut is
explained in item (ii) of Sec. IV.2. The number of lattice points in this
interval is $n_{p}=18$ for the setup L08 and $n_{p}=16$ for the setup L07. The
quantity $\epsilon_{Q_{3}(p_{i}^{2})}$ represents the error in the
corresponding value of $Q_{3}(p_{i}^{2})$. Then, the use of Eq. (48) yields
$\chi^{2}=119$ for the setup L08, and $\chi^{2}=76$ for the setup L07.

At this point, we compute the probability, $P_{Q_{3}^{0}}$, that our result
for $Q_{3}$ is consistent with the null hypothesis through

$P_{Q_{3}^{0}}=\int_{\chi^{2}}^{\infty}\chi^{2}_{{\rm{\scriptscriptstyle
PDF}}}(n_{p},x)dx=\frac{\Gamma(n_{p}/2,\chi^{2}/2)}{\Gamma(n_{p}/2)}\,,$ (49)

where $\chi^{2}_{{\rm{\scriptscriptstyle PDF}}}$ is the $\chi^{2}$ probability
distribution function. Then, using the above quoted values of $n_{p}$ and
$\chi^{2}$ we obtain $P_{Q_{3}^{0}}=6.48\times 10^{-17}$ for L08, and
$P_{Q_{3}^{0}}=8.68\times 10^{-10}$ for L07. In terms of confidence levels to
discard the null hypothesis, these probabilities translate to $8.4\sigma$ and
$6.1\sigma$, respectively.

##  VI Discussion and Conclusions

The non-Abelian implementation of the SM is associated with the appearance of
longitudinally coupled massless poles in the fundamental vertices of the
theory. In addition to endowing the gluons with an effective mass, the action
of the SM leads to the displacement of the Ward identities satisfied by the
fundamental vertices of the theory by amounts proportional to the momentum-
dependent residues of these poles. This property is particularly crucial,
because it enables us to confirm the nonvanishing nature of the pole residues
by forming the difference between the lattice data and the WI prediction for
certain vertex form factors in the soft-gluon limit. In the present work we
have focused on the manifestation of this characteristic effect at the level
of the quark-gluon vertex of QCD with two light quark flavors, $N_{f}=2$.

Our analysis demonstrates that the quark-gluon vertex contains massless poles:
they correspond to the propagators of scalar, color-carrying bound states,
which are composed out of two gluons or of a quark-antiquark pair. The residue
functions associated with these poles play the role of the BS amplitudes for
the bound-state formation, and are determined from a coupled system of BS
equations. Crucially, this system involves also the unquenched version of the
function ${\mathbb{C}}(p^{2})$, known from the well-studied case of the
(quenched) three-gluon vertex.

The residue function $Q_{3}(p^{2})$, associated with the tree-level Dirac
structure $\gamma_{\alpha}$, plays a prominent role in this study, because it
is directly responsible for the displacement produced between the lattice data
and the WI prediction for the classical form factor of the quark-gluon vertex.
In fact, the detailed comparison with the lattice data of Kızılersü  _et al._
(2021) reveals a statistically significant discrepancy (signal), whose
momentum-dependence is completely consistent with the BS solution for
$Q_{3}(p^{2})$, as may be seen in Fig. 10. This finding further corroborates
the action of the SM in the context of unquenched QCD, providing additional
evidence for the robustness of this particular mass generating scenario.

The robustness of the result shown in Fig. 10 hinges on one’s ability to
accurately determine the WI prediction for the form factor
$\lambda_{1}(p^{2})$, namely the $\lambda_{i}^{\star}(p^{2})$ given by Eq.
(27). There, in addition to the components $A(p^{2})$ and $B(p^{2})$ of the
quark propagator, the form factors $K_{1}(p^{2})$ and $K_{4}(p^{2})$,
associated with the quark-ghost kernel [see Eq. (23) and Fig. 4], play a
central role. The determination of $K_{1}(p^{2})$ and $K_{4}(p^{2})$ involves
the fully-dressed quark-gluon vertex, which has been approximated by the
Ansatz given in the second line of Eq. (36). This particular form is inspired
by the special property of “planar degeneracy”, satisfied at a high level of
accuracy by the three-gluon vertex Huber (2020a); Eichmann  _et al._ (2014);
Blum  _et al._ (2014); Pinto-Gómez  _et al._ (2023); Aguilar  _et al._
(2023c). In the case of the quark-gluon vertex, the veracity of this Ansatz
has not been confirmed at a corresponding level of accuracy, even though it
seems fairly plausible, given the observations made below Eq. (36). To be
sure, a more detailed analysis of this entire issue is required, in order to
determine the limitations of this approximation, and the possible
modifications that may arise in the form factors $K_{1}(p^{2})$ and
$K_{4}(p^{2})$. We hope to embark on such a study in the near future.

Finally, even though we have focused on the displacement produced to
$\lambda_{1}(p^{2})$ by $Q_{3}(p^{2})$, the form factors $\lambda_{2}(p^{2})$
and $\lambda_{3}(p^{2})$ get also displaced by the functions
${\mathbb{Q}}_{2+3}(p^{2})$ and ${\mathbb{Q}}_{1}(p^{2})$, respectively [viz.
second and third relation in Eq. (33)]. In fact, a preliminary analysis
indicates that in the case of $\lambda_{3}(p^{2})$ a detectable displacement
may be obtained. This happens despite the fact that ${\mathbb{Q}}_{1}(p^{2})$
is significantly smaller compared to ${\mathbb{C}}(r^{2})$ and $Q_{3}(p^{2})$,
because $\lambda_{3}(p^{2})$ is significantly smaller to $\lambda_{1}(p^{2})$,
and the available lattice data display rather small error bars. It is clearly
important to explore this possibility further, establishing the statistical
significance of the resulting displacement and its similarity to the BS
prediction.

##  VII Acknowledgments

We thank J. Rodríguez-Quintero and O. Oliveira for useful communications. The
work of A. C. A. is supported by the CNPq grant 307854/2019-1. A. C. A also
acknowledges financial support from project 464898/2014-5 (INCT-FNA). D. I.
acknowledges the support provided by the University Centre EDEM. M. N. F. and
J. P. are supported by the Spanish MICINN grant PID2020-113334GB-I00. M. N. F.
acknowledges financial support from Generalitat Valenciana through contract
CIAPOS/2021/74. J. P. also acknowledges funding from the Generalitat
Valenciana grant CIPROM/2022/66.

##  Appendix A Input fits

In this Appendix we present the fits to lattice data for the propagators and
vertex form factors used as external inputs in our calculations. These fits
capture a variety of physical properties in the infrared, discussed elsewhere
Aguilar  _et al._ (2020, 2021a); Ferreira and Papavassiliou (2023), and
reproduce the one-loop anomalous dimensions of the corresponding functions in
the ultraviolet Davydychev  _et al._ (1996, 2001).

The fitting functions employed share certain basic building blocks, which
differ only in the specific values that some sets of parameters acquire in
each case. In particular, the parameters $\delta$, $\kappa$, $\eta_{i}$, and
$b_{i}$ appear in various fits, taking distinct numerical values in each case.
The sets of numerical values employed are displayed in Table 2.

| $\delta$ |  $\kappa^{2}$ [GeV2] |  $\eta_{1}^{2}$ [GeV2] |  $\eta_{2}^{2}$ [GeV2] | $b_{0}$ |  $b_{1}$ [GeV2] |  $b_{2}$ [GeV2] |  $b_{3}$ [GeV2]
---|---|---|---|---|---|---|---|---
$\Delta(r^{2})$ | $0.112$ | $71.8$ | $10.1$ | $0.895$ | $-0.0998$ | $-1.67$ | $0.684$ | $0.321$
$F(r^{2})$ | – | – | $2.68$ | $1.27$ | $-0.440$ | $6.23$ | $0.820$ | $23.2$
$\mathit{L}_{{sg}}(r^{2})$ | $0.0629$ | $12.3$ | $1.00$ | $1.48$ | $0.102$ | $25.9$ | $1.70$ | $19.0$
$A(p^{2})$ [L08] | – | $0.930$ | – | – | $0.360$ | $0.642$ | $0.175$ | $0.462$
$A(p^{2})$ [L07] | – | $2.69$ | – | – | $0.849$ | $0.351$ | $0.193$ | $1.98$
${\cal M}(p^{2})$ [L08] | $0.294$ | $0.520$ | $1.31$ | $13.0$ | – | – | – | –
${\cal M}(p^{2})$ [L07] | $0.110$ | $0.596$ | $1.00$ | $31.1$ | – | – | – | –
$\lambda_{1}(p^{2})$ [L08] | – | – | $24.6$ | $0.0371$ | $-0.587$ | $11.8$ | $0.599$ | $1.03$
$\lambda_{1}(p^{2})$ [L07] | – | – | $24.6$ | $0.0371$ | $-0.483$ | $2.33$ | $0.495$ | $0.745$
Table 2:  Fitting parameters used for Eqs. (50), (53), (54), and (55). In
addition, for ${\cal M}(p^{2})$ we use $m_{q}=6.2$ MeV, and $8$ MeV for the
setups L08 and L07, respectively, and $m_{0}=345$ MeV for both. Lastly, we
employ $\Lambda=610$ MeV for all quantities.

For the unquenched quantities $\Delta(r^{2})$, $F(r^{2})$ and
$\mathit{L}_{{sg}}(r^{2})$, we employ the functional forms

$\displaystyle\Delta^{-1}(r^{2})=$
$\displaystyle\,r^{2}\left[\frac{\delta}{1+(r^{2}/\kappa^{2})}\ln\left(\frac{r^{2}}{\mu^{2}}\right)+U^{d_{{\scriptscriptstyle
A}}}(r^{2})\right]+\nu^{2}R(r^{2})\,,\quad$ $\displaystyle
d_{{\scriptscriptstyle A}}=$ $\displaystyle\,\frac{39-4N_{f}}{6\beta_{0}}\,,$
$\displaystyle F^{-1}(r^{2})=$ $\displaystyle\,U^{d_{{\scriptscriptstyle
c}}}(r^{2})+R(r^{2})\,,\quad$ $\displaystyle d_{{\scriptscriptstyle c}}=$
$\displaystyle\,\frac{9}{4\beta_{0}}\,,$
$\displaystyle\mathit{L}_{{sg}}(r^{2})=$
$\displaystyle\,1.16\left\\{\frac{\delta}{1+(r^{2}/\kappa^{2})}\ln\left(\frac{r^{2}}{\mu^{2}}\right)+U^{d_{{\scriptscriptstyle
3g}}}(r^{2})+R(r^{2})\right\\}\,,\quad$ $\displaystyle d_{{\scriptscriptstyle
3g}}=$ $\displaystyle\,\frac{51-8N_{f}}{12\beta_{0}}\,,$ (50)

where we define $\beta_{0}=11-2N_{f}/3$,

$U(r^{2}):=1+\frac{\ln\left[(r^{2}+\eta^{2}(r^{2}))/(\mu^{2}+\eta^{2}(r^{2}))\right]}{\ln\left(\mu^{2}/\Lambda^{2}\right)}\,,\qquad\eta^{2}(r^{2}):=\frac{\eta_{1}^{2}}{1+r^{2}/\eta_{2}^{2}}\,,$
(51)
$R(r^{2}):=\frac{b_{0}+r^{2}/b_{1}^{2}}{1+r^{2}/b_{2}^{2}+(r^{2}/b_{3}^{2})^{2}}-\frac{b_{0}+\mu^{2}/b_{1}^{2}}{1+\mu^{2}/b_{2}^{2}+(r^{2}/b_{3}^{2})^{2}}\,,$
(52)

and the parameter $\nu=1$ GeV makes the dimensionalities of $R(r^{2})$ and
$\Delta^{-1}(r^{2})$ consistent without changing the dimensions of the
parameters $b_{i}$. Note that the factor of $1.16$ in
$\mathit{L}_{{sg}}(r^{2})$ converts the lattice data of Aguilar  _et al._
(2020) to the $\widetilde{\text{MOM}}$ scheme [see Eq. (72)].

The lattice data of Oliveira  _et al._ (2019b); Kızılersü  _et al._ (2021) for
the quark propagator are given in terms of the quark wave function,
$A^{-1}(p^{2})$, and running mass, ${\cal M}(p^{2}):=B(p^{2})/A(p^{2})$. These
are fitted with

$A(p^{2})=\frac{T(p^{2})}{T(\mu^{2})}\,,\qquad
T(p^{2}):=\frac{b_{0}+p^{2}/b_{1}+(p^{2}/\kappa^{2})^{2}}{1+p^{2}/b_{2}+(p^{2}/b_{3})^{2}}\,,$
(53)

and

${\cal
M}(p^{2})=\frac{m_{0}}{1+(p^{2}/\kappa^{2})^{1+\delta}}+m_{q}\left[\frac{1}{2}\ln\left(\frac{p^{2}+\eta^{2}(p^{2})}{\Lambda^{2}+\eta^{2}(p^{2})}\right)\right]^{-d_{{\scriptscriptstyle{\cal
M}}}}\,,\qquad d_{{\scriptscriptstyle{\cal M}}}=\frac{4}{\beta_{0}}\,.$ (54)

Finally, the fit for $\lambda_{1}(p^{2})$ is obtained with

$\lambda_{1}(p^{2})=\left[U^{d_{{\scriptscriptstyle
qg}}}(p^{2})+R(p^{2})\right]^{-1}\,,\qquad d_{{\scriptscriptstyle
qg}}=\frac{9}{4\beta_{0}}\,.$ (55)

The values of the fitting parameters are collected in Table 2 and its caption.
The resulting curves for $\Delta(r^{2})$ and $\lambda_{1}(p^{2})$ are compared
to the corresponding lattice data in the left panels of Figs. 8 and 10; the
functions $F(r^{2})$, $\mathit{L}_{{sg}}(r^{2})$, $A^{-1}(p^{2})$, and ${\cal
M}(p^{2})$ are shown in Fig. 11.

Figure 11:  Top: lattice data (points) for $F(r^{2})$ from Ayala  _et al._
(2012); Binosi  _et al._ (2017b) (left) and for $\mathit{L}_{{sg}}(r^{2})$
from Aguilar  _et al._ (2020) (right); the corresponding fits (blue
continuous) are given by Eq. (50). Note that for $\mathit{L}_{{sg}}(r^{2})$ we
use $N_{f}=2+1$ data, as discussed in item the (iii) of Sec. IV.2. Bottom:
lattice data of Oliveira  _et al._ (2019b); Kızılersü  _et al._ (2021) for the
quark wave function (left) and running mass (right), and the fits of Eqs. (53)
and (54). The lattice data for the setups L08 and L07 of Oliveira  _et al._
(2019b); Kızılersü  _et al._ (2021) are represented by blue and red points,
respectively, and their fits by blue continuous and red dashed curves.

##  Appendix B Changing renormalization schemes

In this Appendix we show how the lattice inputs given in different
renormalization schemes may be converted to the $\widetilde{\text{MOM}}$,
defined in Eq. (41).

We begin by introducing the renormalization constants

$\displaystyle\Delta_{{\rm{\scriptscriptstyle B}}}$
$\displaystyle=Z_{A}\Delta_{{\rm{\scriptscriptstyle R}}}\,,$ $\displaystyle
F_{\\!{\rm{\scriptscriptstyle B}}}$
$\displaystyle=Z_{c}F_{\\!{\rm{\scriptscriptstyle R}}}\,,$ $\displaystyle
A_{{\rm{\scriptscriptstyle B}}}$ $\displaystyle=Z_{\rm
F}^{-1}A_{{\rm{\scriptscriptstyle R}}}\,,$
$\displaystyle{\rm{I}}\\!\Gamma^{\mu}_{\\!\\!c,\,{\rm{\scriptscriptstyle B}}}$
$\displaystyle=Z_{1}^{-1}{\rm{I}}\\!\Gamma^{\mu}_{\\!\\!c,\,{{\rm{\scriptscriptstyle
R}}}}\,,$ $\displaystyle{\rm{I}}\\!\Gamma_{\\!\\!{{\rm{\scriptscriptstyle
B}}}}^{\mu}$
$\displaystyle=Z_{2}^{-1}{\rm{I}}\\!\Gamma^{\mu}_{\\!\\!{{\rm{\scriptscriptstyle
R}}}}\,,$ $\displaystyle{\rm{I}}\\!\Gamma_{\\!\\!{{\rm{\scriptscriptstyle
B}}}}^{\alpha\mu\nu}$
$\displaystyle=Z_{3}^{-1}{\rm{I}}\\!\Gamma_{\\!\\!{{\rm{\scriptscriptstyle
R}}}}^{\alpha\mu\nu}\,,$ $\displaystyle H_{{\rm{\scriptscriptstyle B}}}$
$\displaystyle=Z_{H}^{-1}H_{{\rm{\scriptscriptstyle R}}}\,,$ $\displaystyle
g_{{\rm{\scriptscriptstyle B}}}$
$\displaystyle=Z_{g}g_{{\rm{\scriptscriptstyle R}}}\,,$ $\displaystyle Z_{g}$
$\displaystyle=Z_{3}Z_{A}^{-3/2}\,,$ (56)

where we denote with subscripts “${\rm B}$” and “${\rm R}$” the bare and
renormalized quantities, respectively. Then, the various STIs of the theory
lead to the relations Marciano and Pagels (1978)

$\frac{Z_{3}}{Z_{A}}=\frac{Z_{1}}{Z_{c}}=\frac{Z_{2}}{Z_{\rm F}}\,.$ (57)

Moreover, in the Landau gauge, we have that $Z_{H}=Z_{1}$, where $Z_{1}$ is
finite Taylor (1971). From now on, the index “R” will be substituted by the
notation introduced below, identifying the different renormalization schemes
under consideration.

Specifically, the required conversions involve three different renormalization
schemes, namely:

  1. 1.

The $\widetilde{\text{MOM}}$ scheme, defined by Eq. (41), which is used in the
lattice determination of $\lambda_{1}(p^{2})$. Quantities without any special
renormalization index, e.g., $\lambda_{1}(p^{2})$, will be understood to be
renormalized in the $\widetilde{\text{MOM}}$ scheme.

  2. 2.

The “asymmetric MOM” scheme,

$\Delta_{{\rm{\scriptscriptstyle asym}}}^{-1}(\mu^{2})=\mu^{2}\,,\quad
F_{\\!\\!{\rm{\scriptscriptstyle asym}}}(\mu^{2})=1\,,\quad
A_{{\rm{\scriptscriptstyle
asym}}}(\mu^{2})=1\,,\quad\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle
asym}}}(\mu^{2})=1\,.$ (58)

which is specified by the index “asym”. This scheme was used in Aguilar  _et
al._ (2020) for the lattice determination of the ($N_{f}=3$) classical form
factor $\mathit{L}_{{sg}}(p^{2})$ of the three-gluon vertex.

  3. 3.

The Taylor scheme, which capitalizes on the Landau gauge finiteness of the
ghost-gluon vertex Taylor (1971), and is defined by the prescription

$\Delta_{{\rm{\scriptscriptstyle T}}}^{-1}(\mu^{2})=\mu^{2}\,,\quad
F_{{\rm{\scriptscriptstyle T}}}(\mu^{2})=1\,,\quad A_{{\rm{\scriptscriptstyle
T}}}(\mu^{2})=1\,,\quad Z^{{\rm{\scriptscriptstyle
T}}}_{1}=Z^{{\rm{\scriptscriptstyle T}}}_{H}=1\,,$ (59)

where we assign an index “${\rm T}$” to quantities renormalized in this
scheme. The value of the coupling constant determined in the Taylor scheme in
Blossier  _et al._ (2010) will be used to determine the coupling constants in
the above two schemes.

Now, let us note that combining Eq. (56) with each of the renormalization
prescriptions of Eqs. (41), (58) and (59) implies

$Z_{A}=Z_{A}^{{\rm{\scriptscriptstyle asym}}}=Z_{A}^{{\rm{\scriptscriptstyle
T}}}=\mu^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})\,,\quad\Longrightarrow\quad\Delta(p^{2})=\Delta_{{\rm{\scriptscriptstyle
asym}}}(p^{2})=\Delta_{{\rm{\scriptscriptstyle T}}}(p^{2})\,,$ (60)

since the bare gluon propagator does not depend on the renormalization scheme.
Similarly, $F(p^{2})=F_{\\!\\!{\rm{\scriptscriptstyle
asym}}}(p^{2})=F_{{\rm{\scriptscriptstyle T}}}(p^{2})$ and
$A(p^{2})=A_{{\rm{\scriptscriptstyle asym}}}(p^{2})=A_{{\rm{\scriptscriptstyle
T}}}(p^{2})$. Hence, we can suppress the indices specifying renormalization
schemes in the propagators.

On the other hand, the renormalization constant of a given vertex has a
different value in each of the above schemes, e.g.,
$Z_{3}^{{\rm{\scriptscriptstyle asym}}}\neq Z_{3}$. As it turns out, the
ratios between these values, e.g., $Z_{3}^{{\rm{\scriptscriptstyle
asym}}}/Z_{3}$, can be conveniently expressed in terms of three special
renormalization group invariant (RGI) effective couplings, defined in terms of
the quark-gluon, three-gluon, and ghost-gluon vertices as

$\displaystyle\alpha_{{\scriptscriptstyle qg}}(p^{2})=$
$\displaystyle\,\alpha_{s}^{{\rm{\scriptscriptstyle
B}}}p^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(p^{2})[\lambda_{1}^{{\rm{\scriptscriptstyle
B}}}(p^{2})]^{2}/A^{2}_{{\rm{\scriptscriptstyle B}}}(p^{2})\,,$
$\displaystyle\alpha_{{\scriptscriptstyle 3g}}(p^{2})=$
$\displaystyle\,\alpha_{s}^{{\rm{\scriptscriptstyle
B}}}[p^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(p^{2})]^{3}[\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle
B}}}(p^{2})]^{2}\,,$ $\displaystyle\alpha_{{\scriptscriptstyle cg}}(p^{2})=$
$\displaystyle\,\alpha_{s}^{{\rm{\scriptscriptstyle
B}}}p^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(p^{2})F^{2}_{{\rm{\scriptscriptstyle B}}}(p^{2})\,.$ (61)

Note that the $\alpha_{{\scriptscriptstyle cg}}(p^{2})$ is the well-known
Taylor coupling Boucaud  _et al._ (2009); Aguilar  _et al._ (2009); Blossier
_et al._ (2010), which takes advantage of the fact that, in Landau gauge, the
unrenormalized ghost-gluon vertex in the soft-ghost configuration reduces to
tree-level Taylor (1971), i.e., $B_{1}^{{\rm{\scriptscriptstyle
B}}}(r^{2},0,r^{2})-B_{2}^{{\rm{\scriptscriptstyle B}}}(r^{2},0,r^{2})=1$, in
the language of Eq. (45).

###  B.1 Taylor to $\widetilde{\text{MOM}}$ conversion and the value of
$Z_{H}$

To convert the value of the coupling between the Taylor and
$\widetilde{\text{MOM}}$ schemes, let us first note that Eqs. (56) and (57)
imply

$Z_{H}=Z_{1}=Z_{c}Z_{2}/Z_{\rm F}\,.$ (62)

Then, using Eq. (41) to write the renormalization constants in terms of the
bare Green’s functions, one shows that

$Z_{H}=\frac{F_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})A_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})}{\lambda_{1}^{{\rm{\scriptscriptstyle
B}}}(\mu^{2})}=\sqrt{\frac{\alpha_{{\scriptscriptstyle
cg}}(\mu^{2})}{\alpha_{{\scriptscriptstyle qg}}(\mu^{2})}}\,.$ (63)

Now, we can use Eq. (63) to determine $\alpha_{{\scriptscriptstyle
qg}}(p^{2})$ from the known $\alpha_{{\scriptscriptstyle cg}}(p^{2})$ of
Blossier  _et al._ (2010). First, at a large momentum scale, here chosen as
$\nu=7$ GeV, we can evaluate Eq. (63) at one loop with massless quarks. Using
the well-known fact that $A_{{\rm{\scriptscriptstyle B}}}(p^{2})=1$ to one
loop in the Landau gauge Pascual and Tarrach (1984); Muta (1987), together
with the corresponding one-loop results for $F_{\\!{\rm{\scriptscriptstyle
B}}}(\nu^{2})$ and $\lambda_{1}^{{\rm{\scriptscriptstyle B}}}(\nu^{2})$ given
in Davydychev  _et al._ (2001), we find

$Z_{H}(\nu=7\leavevmode\nobreak\
\text{GeV})=1+\frac{3\alpha_{s}(\nu=7\leavevmode\nobreak\
\text{GeV})}{16\pi}\,.$ (64)

Then, using the Taylor scheme value of Blossier  _et al._ (2010) for
$\alpha_{s}(\nu=7\leavevmode\nobreak\ \text{GeV})=0.223$, we obtain that
$Z_{H}(\nu=7\leavevmode\nobreak\ \text{GeV})=1.013$, and therefore
$\alpha_{{\scriptscriptstyle qg}}(\nu=7\leavevmode\nobreak\
\text{GeV})=0.217$.

Then, since $\alpha_{{\scriptscriptstyle qg}}(p^{2})$ is RGI, we can rewrite
Eq. (61) in terms of quantities renormalized in the $\widetilde{\text{MOM}}$
scheme with $\mu=2$ GeV, i.e.,

$\alpha_{{\scriptscriptstyle
qg}}(p^{2})=\alpha_{s}p^{2}\Delta(p^{2})\lambda_{1}^{2}(p^{2})/A^{2}(p^{2})\,,$
(65)

and we determine the corresponding value of

$\alpha_{s}(\mu=2\leavevmode\nobreak\ \text{GeV})=0.470(7)\quad[{\rm
L08}],\qquad\alpha_{s}(\mu=2\leavevmode\nobreak\
\text{GeV})=0.471(8)\quad[{\rm L07}],$ (66)

by requiring that Eq. (65) reproduces the previously determined value of
$\alpha_{{\scriptscriptstyle qg}}(\nu=7\leavevmode\nobreak\ \text{GeV})$.

Lastly, with the full momentum dependence of $\alpha_{{\scriptscriptstyle
qg}}(p^{2})$ in hand, we return to Eq. (63) to determine

$Z_{H}=1.120(8)\,\quad[{\rm L08}],\qquad Z_{H}=1.121(9)\,\quad[{\rm L07}].$
(67)

###  B.2 Asymmetric MOM to $\widetilde{\text{MOM}}$ conversion

Now we consider the conversion of $\mathit{L}_{{sg}}(p^{2})$ from the
asymmetric to the $\widetilde{\text{MOM}}$ scheme.

To this end, note first that Eq. (56) implies

$\mathit{L}_{{sg}}(p^{2})=\frac{Z_{3}}{Z_{3}^{{\rm{\scriptscriptstyle
asym}}}}\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle asym}}}(p^{2})\,.$ (68)

Then, our task amounts to computing the ratio
$Z_{3}/Z_{3}^{{\rm{\scriptscriptstyle asym}}}$.

Next, it follows from the renormalization condition of Eq. (58) that
$Z_{3}^{{\rm{\scriptscriptstyle
asym}}}=1/\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle B}}}(\mu^{2})$, whereas
Eqs. (41) and (57) imply

$Z_{3}=\frac{\mu^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})A_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})}{\lambda_{1}^{{\rm{\scriptscriptstyle B}}}(\mu^{2})}\,.$ (69)

Hence, by combining the above results we obtain

$\displaystyle\frac{Z_{3}}{Z_{3}^{{\rm{\scriptscriptstyle asym}}}}=$
$\displaystyle\,\frac{\mu^{2}\Delta_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})A_{{\rm{\scriptscriptstyle
B}}}(\mu^{2})\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle
B}}}(\mu^{2})}{\lambda_{1}^{{\rm{\scriptscriptstyle
B}}}(\mu^{2})}=\sqrt{\frac{\alpha_{{\scriptscriptstyle
3g}}(\mu^{2})}{\alpha_{{\scriptscriptstyle qg}}(\mu^{2})}}\,.$ (70)

Now, to fix the value of $\alpha_{{\scriptscriptstyle 3g}}(\mu^{2})$ we use
the same strategy used in the determination of $\alpha_{{\scriptscriptstyle
qg}}(p^{2})$.

At a large momentum scale, $\nu=7$ GeV, we can evaluate the combination of
bare Green’s functions at one loop. Using the textbook result for
$\Delta_{{\rm{\scriptscriptstyle B}}}(\nu^{2})$ and
$A_{{\rm{\scriptscriptstyle B}}}(\nu^{2})=1$ Pascual and Tarrach (1984); Muta
(1987) and the results of Davydychev  _et al._ (1996, 2001) for the vertices
we obtain

$\left.\frac{Z_{3}}{Z_{3}^{{\rm{\scriptscriptstyle
asym}}}}\right|_{\nu=7\leavevmode\nobreak\
\text{GeV}}\\!\\!\\!=1+\frac{\alpha_{s}(\nu=7\leavevmode\nobreak\
\text{GeV})}{96\pi}(129-16N_{f})=1.072\,,$ (71)

which implies $\alpha_{{\scriptscriptstyle 3g}}(\nu=7\leavevmode\nobreak\
\text{GeV})=0.249$.

Then, from the RGI nature of $\alpha_{{\scriptscriptstyle 3g}}(p^{2})$, we can
compute its entire momentum dependence using ingredients renormalized in the
asymmetric scheme. In particular, we obtain the three-gluon coupling
$\alpha_{{\scriptscriptstyle 3g}}(\mu=2\leavevmode\nobreak\
\text{GeV})=0.635$, and

$\frac{Z_{3}}{Z_{3}^{{\rm{\scriptscriptstyle asym}}}}=1.16\,.$ (72)

##  Appendix C Unquenched ghost-gluon vertex

In this Appendix we briefly describe the determination of the classical form
factor $B_{1}(r^{2},k^{2},q^{2})$ of the ghost-gluon vertex [see Eq. (45)],
which is relevant for the calculation of the quark-ghost form factors,
$K_{i}$.

The one-loop dressed SDE of the ghost-gluon vertex is formally identical to
its quenched version, given diagrammatically in Fig. 12 of Ferreira and
Papavassiliou (2023); the dependence on $N_{f}$ enters through the various
quantities comprising the diagrams, such as the coupling, the gluon and ghost
propagators, and the three-gluon vertex.

For the three-gluon vertex, we use the approximation given in Eq. (36). Then,
an appropriate tensor projection of the SDE yields a dynamical equation for
$B_{1}(r^{2},k^{2},q^{2})$, which is given by Eq. (88) of Ferreira and
Papavassiliou (2023). Finally, this equation is solved numerically, using Eq.
(66) for $\alpha_{s}$, the L08 setup value given in Eq. (67) for
$Z_{1}=Z_{H}$, and the fits given by Eq. (50) for $\Delta(r^{2})$, $F(r^{2})$
and $\mathit{L}_{{sg}}(r^{2})$.

To expedite the comparison between $N_{f}=2$ and quenched results for the
ghost-gluon vertex, we will show its classical form factor renormalized in the
Taylor scheme, $B_{1}^{{\rm{\scriptscriptstyle T}}}$. Its
$\widetilde{\text{MOM}}$ value, $B_{1}$, is obtained immediately through

$B_{1}(r^{2},k^{2},q^{2})=Z_{H}B_{1}^{{\rm{\scriptscriptstyle
T}}}(r^{2},k^{2},q^{2})\,,$ (73)

where $Z_{H}$ is given by Eq. (67).

Figure 12:  Left: Classical form factor, $B_{1}^{{\rm{\scriptscriptstyle
T}}}(r^{2},k^{2},q^{2})$, of the ghost-gluon vertex renormalized in the Taylor
scheme, for a fixed value $\phi=2\pi/3$ of the angle between $k$ and $q$. The
orange continuous line highlights the soft-ghost ($k=0$) limit, respectively.
Right: Comparison between the $N_{f}=2$ (orange continuous) and quenched (blue
dashed) values of the soft-ghost limit, $B_{1}^{{\rm{\scriptscriptstyle
T}}}(q^{2},0,q^{2})$.

On the left panel of Fig. 12 we show $B_{1}^{{\rm{\scriptscriptstyle
T}}}(r^{2},k^{2},q^{2})$ as a function of the magnitudes of the ghost and
gluon momenta, $k$ and $q$, respectively; for the angle $\phi$ between them we
choose the representative value $\phi=2\pi/3$. On this panel we note that
$B_{1}^{{\rm{\scriptscriptstyle T}}}(r^{2},k^{2},q^{2})$ has the same
qualitative behavior as its quenched form obtained in various studies555For
related studies within the Curci-Ferrari model, see Mintz  _et al._ (2018);
Barrios  _et al._ (2020). Schleifenbaum  _et al._ (2005); Ilgenfritz  _et al._
(2007); Aguilar  _et al._ (2013); Cyrol  _et al._ (2016); Huber (2020b);
Aguilar  _et al._ (2019); Huber (2020a); Aguilar  _et al._ (2021b). In
particular, it is nearly identical in shape to the result shown in Fig. 13 of
Ferreira and Papavassiliou (2023), which uses the same approximation of Eq.
(36) for the three-gluon vertex.

Next, we focus on the soft-ghost limit ($k=0$), which appears as an ingredient
in Eq. (46) for the $K_{i}(p^{2})$. This limit corresponds to the slice
highlighted as an orange continuous curve on the left panel of Fig. 12. On the
right panel of the same figure, we compare the $N_{f}=2$ value of
$B_{1}^{{\rm{\scriptscriptstyle T}}}(q^{2},0,q^{2})$ to its quenched
counterpart computed in Ferreira and Papavassiliou (2023), represented here as
a blue dashed line. The results are found to differ by less than $2.75\%$
within the entire momentum range, and are virtually indistinguishable for
$q<1$ GeV. Their most noticeable difference is in the ultraviolet, where the
$N_{f}=2$ result is seen to be systematically larger than its quenched
version. The same pattern is found for all kinematic configurations.

Note that the observed enhancement of the unquenched
$B_{1}^{{\rm{\scriptscriptstyle T}}}(q^{2},0,q^{2})$ in the ultraviolet is
compatible with perturbation theory. Indeed, a one-loop calculation yields

$B_{1}^{{\rm{\scriptscriptstyle
T}}}(q^{2},0,q^{2})=1+\frac{33\alpha_{s}^{{\rm{\scriptscriptstyle
T}}}}{32\pi}\,,$ (74)

which depends on $N_{f}$ only through the value of
$\alpha_{s}^{{\rm{\scriptscriptstyle T}}}$. Then, since
$\alpha_{s}^{{\rm{\scriptscriptstyle T}}}$ increases with $N_{f}$ Ayala  _et
al._ (2012), so does $B_{1}^{{\rm{\scriptscriptstyle T}}}(q^{2},0,q^{2})$.

## References

  * Marciano and Pagels (1978) W. J. Marciano and H. Pagels, Phys. Rept. 36, 137 (1978).
  * Schwinger (1962a) J. S. Schwinger, Phys. Rev. 125, 397 (1962a).
  * Schwinger (1962b) J. S. Schwinger, Phys. Rev. 128, 2425 (1962b).
  * Smit (1974) J. Smit, Phys. Rev. D 10, 2473 (1974).
  * Eichten and Feinberg (1974) E. Eichten and F. Feinberg, Phys. Rev. D 10, 3254 (1974).
  * Cornwall (1982) J. M. Cornwall, Phys. Rev. D 26, 1453 (1982).
  * Aguilar  _et al._ (2012a) A. C. Aguilar, D. Ibanez, V. Mathieu, and J. Papavassiliou, Phys. Rev. D85, 014018 (2012a).
  * Binosi  _et al._ (2012) D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D86, 085033 (2012).
  * Aguilar  _et al._ (2016a) A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D94, 045002 (2016a).
  * Eichmann  _et al._ (2021) G. Eichmann, J. M. Pawlowski, and J. a. M. Silva, Phys. Rev. D 104, 114016 (2021).
  * Aguilar  _et al._ (2022) A. C. Aguilar, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 105, 014030 (2022).
  * Aguilar  _et al._ (2023a) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, F. Pinto-Gómez, C. D. Roberts, and J. Rodríguez-Quintero, Phys. Lett. B 841, 137906 (2023a).
  * Papavassiliou (2022) J. Papavassiliou, Chin. Phys. C 46, 112001 (2022).
  * Ferreira and Papavassiliou (2023) M. N. Ferreira and J. Papavassiliou, Particles 6, 312 (2023).
  * Fischer  _et al._ (2009) C. S. Fischer, A. Maas, and J. M. Pawlowski, Annals Phys. 324, 2408 (2009).
  * Wilson  _et al._ (1994) K. G. Wilson, T. S. Walhout, A. Harindranath, W.-M. Zhang, R. J. Perry, and S. D. Glazek, Phys. Rev. D49, 6720 (1994).
  * Kondo (2001) K.-I. Kondo, Phys. Lett. B 514, 335 (2001).
  * Epple  _et al._ (2008) D. Epple, H. Reinhardt, W. Schleifenbaum, and A. P. Szczepaniak, Phys. Rev. D77, 085007 (2008).
  * Campagnari and Reinhardt (2010) D. R. Campagnari and H. Reinhardt, Phys. Rev. D82, 105021 (2010).
  * Tissier and Wschebor (2010) M. Tissier and N. Wschebor, Phys. Rev. D 82, 101701 (2010).
  * Alkofer and Alkofer (2011) N. Alkofer and R. Alkofer, Phys. Lett. B 702, 158 (2011).
  * Serreau and Tissier (2012) J. Serreau and M. Tissier, Phys. Lett. B712, 97 (2012).
  * Siringo (2016) F. Siringo, Nucl. Phys. B907, 572 (2016).
  * Cyrol  _et al._ (2016) A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D94, 054005 (2016).
  * Glazek  _et al._ (2017) S. D. Glazek, M. Gómez-Rocha, J. More, and K. Serafin, Phys. Lett. B773, 172 (2017).
  * Huber (2020a) M. Q. Huber, Phys. Rev. D 101, 114009 (2020a).
  * Horak  _et al._ (2022) J. Horak, F. Ihssen, J. Papavassiliou, J. M. Pawlowski, A. Weber, and C. Wetterich, SciPost Phys. 13, 042 (2022).
  * Zumino (1965) B. Zumino, Acta Phys. Austriaca Suppl. 2, 212 (1965).
  * Jackiw and Johnson (1973) R. Jackiw and K. Johnson, Phys. Rev. D 8, 2386 (1973).
  * Jackiw (1973) R. Jackiw, In *Erice 1973, Proceedings, Laws Of Hadronic Matter*, New York 1975, 225-251 and M I T Cambridge - COO-3069-190 (73,REC.AUG 74) 23p (1973).
  * Cornwall and Norton (1973) J. Cornwall and R. Norton, Phys. Rev. D 8, 3338 (1973).
  * Poggio  _et al._ (1975) E. C. Poggio, E. Tomboulis, and S. H. H. Tye, Phys. Rev. D11, 2839 (1975).
  * Cornwall (1979) J. M. Cornwall, Nucl. Phys. B 157, 392 (1979).
  * Aguilar  _et al._ (2008) A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D78, 025010 (2008).
  * Binosi and Papavassiliou (2018) D. Binosi and J. Papavassiliou, Phys. Rev. D97, 054029 (2018).
  * Aguilar  _et al._ (2018) A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Eur. Phys. J. C78, 181 (2018).
  * Roberts and Williams (1994) C. D. Roberts and A. G. Williams, Prog. Part. Nucl. Phys. 33, 477 (1994).
  * Cloet and Roberts (2014) I. C. Cloet and C. D. Roberts, Prog. Part. Nucl. Phys. 77, 1 (2014).
  * Alkofer and von Smekal (2001) R. Alkofer and L. von Smekal, Phys. Rept. 353, 281 (2001).
  * Fischer (2006) C. S. Fischer, J. Phys. G 32, R253 (2006).
  * Roberts (2008) C. D. Roberts, Prog. Part. Nucl. Phys. 61, 50 (2008).
  * Binosi and Papavassiliou (2009) D. Binosi and J. Papavassiliou, Phys. Rept. 479, 1 (2009).
  * Aguilar  _et al._ (2016b) A. C. Aguilar, D. Binosi, and J. Papavassiliou, Front. Phys.(Beijing) 11, 111203 (2016b).
  * Huber (2020b) M. Q. Huber, Phys. Rept. 879, 1 (2020b).
  * Athenodorou  _et al._ (2016) A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, Phys. Lett. B761, 444 (2016).
  * Duarte  _et al._ (2016) A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D94, 074502 (2016).
  * Boucaud  _et al._ (2017) P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D95, 114503 (2017).
  * 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).
  * 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. C80, 154 (2020).
  * Aguilar  _et al._ (2021a) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, and J. Rodríguez-Quintero, Phys. Lett. B 818, 136352 (2021a).
  * Aguilar  _et al._ (2021b) A. C. Aguilar, C. O. Ambrósio, F. De Soto, M. N. Ferreira, B. M. Oliveira, J. Papavassiliou, and J. Rodríguez-Quintero, Phys. Rev. D 104, 054028 (2021b).
  * Ball and Chiu (1980) J. S. Ball and T.-W. Chiu, Phys. Rev. D22, 2542 (1980).
  * Kizilersu  _et al._ (1995) A. Kizilersu, M. Reenders, and M. Pennington, Phys. Rev. D52, 1242 (1995).
  * Davydychev  _et al._ (2001) A. I. Davydychev, P. Osland, and L. Saks, Phys. Rev. D63, 014022 (2001).
  * Gracey (2014) J. Gracey, Phys. Rev. D 90, 025014 (2014).
  * Gracey (2011) J. A. Gracey, Phys. Rev. D84, 085011 (2011).
  * Bermudez  _et al._ (2017) R. Bermudez, L. Albino, L. X. Gutiérrez-Guerrero, M. E. Tejeda-Yeomans, and A. Bashir, Phys. Rev. D95, 034041 (2017).
  * Bender  _et al._ (1996) A. Bender, C. D. Roberts, and L. Von Smekal, Phys. Lett. B380, 7 (1996).
  * Bhagwat and Tandy (2004) M. Bhagwat and P. Tandy, Phys. Rev. D70, 094039 (2004).
  * Llanes-Estrada  _et al._ (2006) F. J. Llanes-Estrada, C. S. Fischer, and R. Alkofer, Nucl. Phys. Proc. Suppl. 152, 43 (2006).
  * Holl  _et al._ (2005) A. Holl, A. Krassnigg, and C. D. Roberts, Nucl. Phys. Proc. Suppl. 141, 47 (2005).
  * Matevosyan  _et al._ (2007) H. H. Matevosyan, A. W. Thomas, and P. C. Tandy, Phys. Rev. C75, 045201 (2007).
  * Aguilar and Papavassiliou (2011) A. C. Aguilar and J. Papavassiliou, Phys. Rev. D83, 014013 (2011).
  * Qin  _et al._ (2013) S.-X. Qin, L. Chang, Y.-X. Liu, C. D. Roberts, and S. M. Schmidt, Phys. Lett. B722, 384 (2013).
  * Binosi  _et al._ (2017a) D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D95, 031501 (2017a).
  * Hopfer  _et al._ (2012) M. Hopfer, A. Windisch, and R. Alkofer, PoS ConfinementX, 073 (2012).
  * Aguilar  _et al._ (2012b) A. C. Aguilar, D. Binosi, J. C. Cardona, and J. Papavassiliou, PoS ConfinementX, 103 (2012b).
  * Rojas  _et al._ (2013) E. Rojas, J. de Melo, B. El-Bennich, O. Oliveira, and T. Frederico, J. High Energy Phys. 10, 193 (2013).
  * Williams (2015) R. Williams, Eur. Phys. J. A51, 57 (2015).
  * Williams  _et al._ (2016) R. Williams, C. S. Fischer, and W. Heupel, Phys. Rev. D93, 034026 (2016).
  * Sanchis-Alepuz and Williams (2015) H. Sanchis-Alepuz and R. Williams, Phys. Lett. B749, 592 (2015).
  * Peláez  _et al._ (2015) M. Peláez, M. Tissier, and N. Wschebor, Phys. Rev. D 92, 045012 (2015).
  * Alkofer  _et al._ (2009) R. Alkofer, C. S. Fischer, F. J. Llanes-Estrada, and K. Schwenzer, Annals Phys. 324, 106 (2009).
  * Mitter  _et al._ (2015) M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D91, 054035 (2015).
  * Cyrol  _et al._ (2018) A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D97, 054006 (2018).
  * Aguilar  _et al._ (2014) A. C. Aguilar, D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D90, 065027 (2014).
  * Gao  _et al._ (2021) F. Gao, J. Papavassiliou, and J. M. Pawlowski, Phys. Rev. D 103, 094013 (2021).
  * Aguilar  _et al._ (2017) A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D96, 014029 (2017).
  * Oliveira  _et al._ (2019a) O. Oliveira, W. de Paula, T. Frederico, and J. P. B. C. de Melo, Eur. Phys. J. C 79, 116 (2019a).
  * Albino  _et al._ (2019) L. Albino, A. Bashir, L. X. G. Guerrero, B. E. Bennich, and E. Rojas, Phys. Rev. D 100, 054028 (2019).
  * Tang  _et al._ (2019) C. Tang, F. Gao, and Y.-X. Liu, Phys. Rev. D 100, 056001 (2019).
  * Albino  _et al._ (2021) L. Albino, A. Bashir, B. El-Bennich, E. Rojas, F. E. Serna, and R. C. da Silveira, JHEP 11, 196 (2021).
  * Skullerud  _et al._ (2002) J. Skullerud, P. O. Bowman, and A. Kizilersu, in  _5th International Conference on Quark Confinement and the Hadron Spectrum_ (2002) pp. 270–272.
  * Skullerud and Kizilersu (2002) J. Skullerud and A. Kizilersu, JHEP 09, 013 (2002).
  * Skullerud  _et al._ (2003) J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, J. High Energy Phys. 04, 047 (2003).
  * Skullerud  _et al._ (2005) J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, Nucl. Phys. Proc. Suppl. 141, 244 (2005).
  * Lin (2006) H.-W. Lin, Phys. Rev. D73, 094511 (2006).
  * 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).
  * Kizilersu  _et al._ (2007) A. Kizilersu, D. B. Leinweber, J.-I. Skullerud, and A. G. Williams, Eur. Phys. J. C50, 871 (2007).
  * Sternbeck  _et al._ (2017) A. Sternbeck, P.-H. Balduf, A. Kizilersu, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. G. Williams, PoS LATTICE2016, 349 (2017).
  * Skullerud  _et al._ (2022) J.-I. Skullerud, A. Kızılersü, O. Oliveira, P. Silva, and A. Sternbeck, PoS LATTICE2021, 305 (2022).
  * Oliveira  _et al._ (2016) O. Oliveira, A. Kizilersu, P. J. Silva, J.-I. Skullerud, A. Sternbeck, and A. G. Williams, Acta Phys. Polon. Supp. 9, 363 (2016).
  * Oliveira  _et al._ (2018) O. Oliveira, T. Frederico, W. de Paula, and J. P. B. C. de Melo, Eur. Phys. J. C78, 553 (2018).
  * Ayala  _et al._ (2012) A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, and J. Rodriguez-Quintero, Phys. Rev. D86, 074512 (2012).
  * Binosi  _et al._ (2017b) D. Binosi, C. D. Roberts, and J. Rodriguez-Quintero, Phys. Rev. D 95, 114009 (2017b).
  * Aguilar  _et al._ (2023b) A. C. Aguilar, M. N. Ferreira, B. M. Oliveira, J. Papavassiliou, and L. R. Santos, [arXiv:2306.16283 [hep-ph]] .
  * Ibañez and Papavassiliou (2013) D. Ibañez and J. Papavassiliou, Phys. Rev. D87, 034008 (2013).
  * Taylor (1971) J. Taylor, Nucl. Phys. B 33, 436 (1971).
  * Slavnov (1972) A. Slavnov, Theor. Math. Phys. 10, 99 (1972).
  * Natale and Rodrigues da Silva (1997) A. A. Natale and P. S. Rodrigues da Silva, Phys. Lett. B 392, 444 (1997).
  * Fleischer  _et al._ (1999) J. Fleischer, F. Jegerlehner, O. V. Tarasov, and O. L. Veretin, Nucl. Phys. B 539, 671 (1999), [Erratum: Nucl.Phys.B 571, 511–512 (2000)].
  * Franco and Lubicz (1998) E. Franco and V. Lubicz, Nucl. Phys. B 531, 641 (1998).
  * Maris and Roberts (2003) P. Maris and C. D. Roberts, Int. J. Mod. Phys. E12, 297 (2003).
  * Bhagwat  _et al._ (2003) M. S. Bhagwat, M. A. Pichowsky, C. D. Roberts, and P. C. Tandy, Phys. Rev. C 68, 015203 (2003).
  * Bowman  _et al._ (2005) P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, A. G. Williams, and J.-b. Zhang, Phys. Rev. D71, 054507 (2005).
  * Furui and Nakajima (2006) S. Furui and H. Nakajima, Phys. Rev. D 73, 074503 (2006).
  * Sauli  _et al._ (2007) V. Sauli, J. Adam, Jr., and P. Bicudo, Phys. Rev. D75, 087701 (2007).
  * Chang  _et al._ (2007) L. Chang, Y.-X. Liu, M. S. Bhagwat, C. D. Roberts, and S. V. Wright, Phys. Rev. C 75, 015201 (2007).
  * Bashir  _et al._ (2012) A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y.-X. Liu, _et al._ , Commun. Theor. Phys. 58, 79 (2012).
  * Oliveira  _et al._ (2019b) O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Phys. Rev. D 99, 094506 (2019b).
  * Chen  _et al._ (2021) L.-f. Chen, Z. Bai, F. Gao, and Y.-x. Liu, Phys. Rev. D 104, 094041 (2021).
  * Eichmann  _et al._ (2014) G. Eichmann, R. Williams, R. Alkofer, and M. Vujinovic, Phys. Rev. D89, 105014 (2014).
  * Blum  _et al._ (2014) A. Blum, M. Q. Huber, M. Mitter, and L. von Smekal, Phys. Rev. D89, 061703 (2014).
  * Aguilar  _et al._ (2023c) A. C. Aguilar, M. N. Ferreira, J. Papavassiliou, and L. R. Santos, Eur. Phys. J. C 83, 549 (2023c).
  * Blum (2017) A. L. Blum, Ph.D. thesis, Graz U. (2017).
  * Blum  _et al._ (2015) A. L. Blum, R. Alkofer, M. Q. Huber, and A. Windisch, Acta Phys. Polon. Supp. 8, 321 (2015).
  * Skullerud  _et al._ (2001) J. Skullerud, D. B. Leinweber, and A. G. Williams, Phys. Rev. D 64, 074508 (2001).
  * Pascual and Tarrach (1984) P. Pascual and R. Tarrach, Lect. Notes Phys. 194, 1 (1984).
  * Muta (1987) T. Muta, World Sci. Lect. Notes Phys. 5, 1 (1987).
  * Davydychev  _et al._ (1996) A. I. Davydychev, P. Osland, and O. Tarasov, Phys. Rev. D 54, 4087 (1996), [Erratum: Phys.Rev.D 59, 109901 (1999)].
  * Blossier  _et al._ (2010) B. Blossier, P. Boucaud, F. De soto, V. Morenas, M. Gravina, O. Pene, and J. Rodriguez-Quintero (ETM), Phys. Rev. D 82, 034510 (2010).
  * Boucaud  _et al._ (2009) P. Boucaud, F. De Soto, J. Leroy, A. Le Yaouanc, J. Micheli, _et al._ , Phys. Rev. D79, 014508 (2009).
  * Aguilar  _et al._ (2009) A. C. Aguilar, D. Binosi, J. Papavassiliou, and J. Rodriguez-Quintero, Phys. Rev. D80, 085018 (2009).
  * Mintz  _et al._ (2018) B. W. Mintz, L. F. Palhares, S. P. Sorella, and A. D. Pereira, Phys. Rev. D97, 034020 (2018).
  * Barrios  _et al._ (2020) N. Barrios, M. Peláez, U. Reinosa, and N. Wschebor, Phys. Rev. D 102, 114016 (2020).
  * Schleifenbaum  _et al._ (2005) W. Schleifenbaum, A. Maas, J. Wambach, and R. Alkofer, Phys. Rev. D 72, 014017 (2005).
  * Ilgenfritz  _et al._ (2007) E.-M. Ilgenfritz, M. Muller-Preussker, A. Sternbeck, A. Schiller, and I. Bogolubsky, Braz.J. Phys. 37, 193 (2007).
  * Aguilar  _et al._ (2013) A. C. Aguilar, D. Ibañez, and J. Papavassiliou, Phys. Rev. D87, 114020 (2013).
  * Aguilar  _et al._ (2019) A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D99, 034026 (2019).