# Nonlinear Schwinger mechanism in QCD,
and Fredholm alternatives theorem

M. N. Ferreira  School of Physics, Nanjing University, Nanjing, Jiangsu
210093, China  Institute for Nonperturbative Physics, Nanjing University,
Nanjing, Jiangsu 210093, China  J. Papavassiliou  Department of Theoretical
Physics and IFIC, University of Valencia and CSIC, E-46100, Valencia, Spain
ExtreMe Matter Institute EMMI, GSI, Planckstrasse 1, 64291, Darmstadt, Germany

###### Abstract

We present a novel implementation of the Schwinger mechanism in QCD, which
fixes uniquely the scale of the effective gluon mass and streamlines
considerably the procedure of multiplicative renormalization. The key
advantage of this method stems from the nonlinear nature of the dynamical
equation that generates massless poles in the longitudinal sector of the
three-gluon vertex. An exceptional feature of this approach is an extensive
cancellation involving the components of the integral expression that
determines the gluon mass; it is triggered once the Schwinger-Dyson equation
of the pole-free part of the three-gluon vertex has been appropriately
exploited. It turns out that this cancellation is driven by the so-called
Fredholm alternatives theorem, which operates among the set of integral
equations describing this system. Quite remarkably, in the linearized
approximation this theorem enforces the exact masslessness of the gluon.
Instead, the nonlinearity induced by the full treatment of the relevant kernel
evades this theorem, allowing for the emergence of a nonvanishing mass. The
numerical results obtained from the resulting equations are compatible with
the lattice findings, and may be further refined through the inclusion of the
remaining fundamental vertices of the theory.

##  I Introduction

As is well-known, the gauge symmetry of the classical Yang-Mills action [1, 2,
3, 4] prohibits the inclusion of a mass term $m^{2}A^{a}_{\mu}A^{a\mu}$ for
the gauge field $A_{\mu}^{a}$, and this prohibition persists after the gauge-
fixing through the Faddeev-Popov procedure [5] has been carried out, due to
the residual Becchi-Rouet-Stora-Tyutin (BRST) symmetry [6, 7, 8]. Moreover,
symmetry-preserving regularization schemes, such as dimensional regularization
[9, 10], prevent the generation of a mass term at any finite order in
perturbation theory. Nevertheless, as has been argued over four decades ago
[11], a non-perturbative effective gluon mass may be generated through the
non-Abelian realization of the Schwinger mechanism (SM) [12, 13]. This
dynamical scenario has received particular attention over the years, becoming
an integral part of the intense activity associated with the emergence of a
mass gap in the gauge sector of QCD [14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53].

In its original formulation, the SM is based on the fundamental observation
that a mass for a vector meson may be generated, without clashing with the
gauge symmetry of the action, if the corresponding vacuum polarization
develops a pole with positive residue ($m^{2}$) at zero-momentum transfer [12,
13]. In that case, the gluon propagator, $\Delta(q)$, saturates to a fixed
value at the origin, namely $\Delta^{-1}(0)=m^{2}$, such that the residue of
the pole plays the role of an effective mass.

The implementation of this basic idea in the context of Yang-Mills theories,
in general, and in the gauge sector of QCD in particular, is technically
rather complicated. Specifically, the required pole is provided by the
propagator $\delta^{ab}/q^{2}$ of a massless color-carrying scalar bound
state, formed out of the fusion of a pair of gluons [14, 15, 54, 11, 55, 56,
57, 58, 59]. The amplitude for the formation of such a pole, denoted by
$\mathbb{B}(r)$, is governed by a homogeneous Bethe-Salpeter equation (BSE).
The existence of a nontrivial solution for $\mathbb{B}(r)$ endows the three-
gluon vertex with a longitudinally coupled pole [55, 56, 57, 58, 59], which is
transmitted to the Schwinger-Dyson equation (SDE) [60, 61, 62, 63, 64, 65, 66,
67, 68, 69, 70, 57, 71, 72, 73, 74, 51, 53] that governs the momentum
evolution of the gluon propagator, triggering the SM. In the limit $q\to 0$,
this SDE expresses the effective gluon mass as an integral involving
$\mathbb{B}(r)$. In addition, the Ward identity (WI) satisfied by the three-
gluon vertex is modified by an amount controlled by the so-called
“displacement function” [75, 59, 76], denoted by ${\mathbb{C}}(r)$, which is
intimately connected to $\mathbb{B}(r)$.

Even though considerable advances have been made in our understanding of the
dynamical scenario outlined above, the treatments presented in the literature
leave certain key questions unanswered. First, in the linearized version of
the BSE studied thus far, the scale of the solutions for $\mathbb{B}(r)$
remains undetermined [55, 56, 57, 58, 59]; thus, the value of the gluon mass
is adjusted from the saturation point of the gluon propagator obtained in
lattice simulations [77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
91, 92, 93, 94]. Second, the integral expression furnishing the mass is
multiplied by $Z_{3}$, the renormalization constant of the three-gluon vertex.
As a result, the full treatment of the gluon mass equation requires the
implementation of multiplicative renormalization at a nonperturbative level;
however, the numerical evaluation of this integral has been carried out by
simply setting $Z_{3}=1$ [55, 40, 56, 57, 59, 95].

In the present work we report significant progress towards the self-consistent
resolution of the open questions mentioned above. The main points of this
analysis may be summarized as follows.

(i) The BSE governing $\mathbb{B}(r)$ is made nonlinear, by including a term
cubic in $\mathbb{B}$, not considered in previous studies [55, 56, 57, 58, 59,
95]. Quite interestingly, the quadratic part of this dependence on
$\mathbb{B}$ is absorbed into a constant denoted by $\omega$; as a result, the
total contribution of this term to the BSE is a factor of the form
$\omega\,\mathbb{B}(r)$. Consequently, the resulting integral equation may
still be treated as an eigenvalue problem, exactly as was done in the past
with the linear BSE. The presence of this novel term determines the scale of
the solutions up to an overall sign; since the basic quantities of interest,
such as the gluon mass and the displacement function, depend quadratically on
$\mathbb{B}$, this residual sign ambiguity is immaterial.

(ii) The multiplicative renormalization is carried out by introducing the SDE
for the three-gluon vertex, which depends on the renormalization constant
$Z_{3}$. It turns out that, by virtue of a massive cancellation triggered by
the form of the BSE in (i), the multiplicative renormalization of the mass
equation may be carried out exactly. The resulting manifestly finite
expression is proportional to the parameter $\omega$. Thus, rather remarkably,
in the absence of the non-cubic term in the BSE, the finite gluon mass
vanishes, even in the presence of a non-vanishing $\mathbb{B}(r)$. Instead,
when $\omega$ assumes its natural value dictated by the BSE, the resulting
gluon mass is non-vanishing; its value can be made compatible with the
$\Delta^{-1}(0)$ obtained from lattice simulations by adjusting the kernel of
the BSE.

(iii) It turns out that the extensive cancellations mentioned in (ii) can be
attributed to a concrete mathematical reason, namely the Fredholm alternatives
theorem [96, 97], which is in full operation when $\omega=0$. In that case,
the kernels of the homogeneous BSE for $\mathbb{B}(r)$ and of the
inhomogeneous SDE for the three-gluon vertex are identical. Then, according to
the theorem, no simultaneous solutions exist to these two integral equations
unless a certain integral expression vanishes; and, quite notably, this
expression is none other than the mass equation derived from the gluon SDE.
Instead, the presence of a non-vanishing $\omega$ introduces a difference
between the two kernels, thus evading the main assumption of the theorem and
leading to the emergence of a gluon mass.

The article is organized as follows. In Sec. II we review the theoretical
framework employed in this work, focusing on the set of special relations
pertaining to the SM. Then, certain key results are derived in Sec. III, which
are extensively used in the ensuing analysis. In Sec. IV the renormalization
of the relevant quantities is outlined, placing particular emphasis on the
treatment of the components emerging from the activation of the SM. In Sec. V
we derive in detail the BSE that controls the formation of the bound state
pole, paying special attention to the novel term that causes the non-linearity
of the resulting integral equation. In Sec. VI we explain in detail how the
BSE of the previous section leads to the fixing of the scale of the solutions
obtained. Then, in Sec. VII we carry out the multiplicative renormalization of
the gluon mass, which proceeds through the implementation of a crucial
cancellation. In Sec. VIII we state the Fredholm alternatives theorem, and
apply it to the set of integral equations considered in the previous section.
In Sec. IX we perform a detailed numerical analysis of the above equations,
and study the dependence of the gluon mass and the displacement function on
the details of the BSE kernel. In Sec. X we discuss our results and summarize
our conclusions. Finally, in Appendix A we illustrate the central mathematical
construction of Sec. VII by means of a concrete example.

##  II General theoretical framework

In this section we introduce the necessary notation, and present a brief
overview of the main notions associated with the implementation of the SM in
the context of Yang-Mills theories; for recent reviews on the subject, see
[51, 53].

Let us consider the (Landau-gauge) gluon propagator

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

and the associated dimensionless vacuum polarization ${\bf\Pi}(q)$,

$\Delta^{-1}({q})=q^{2}[1+i{\bf\Pi}(q)]\,.$ (2)

In addition, the three-gluon vertex, ${\cal G}^{abc}_{\alpha\mu\nu}(q,r,p)$,
with all momenta incoming and $q+r+p=0$, is defined as

${\cal G}^{abc}_{\alpha\mu\nu}(q,r,p)\ =\ \langle
0|T[\widetilde{A}^{a}_{\alpha}(q)\widetilde{A}^{b}_{\mu}(r)\widetilde{A}^{c}_{\nu}(p)]|0\rangle\,,$
(3)

where $\widetilde{A}^{a}_{\alpha}$ denotes the Fourier transformed SU(3) gauge
field, and $T$ the standard time-ordering operation. It is convenient for the
analysis that follows to explicitly factor out of ${\cal
G}^{abc}_{\alpha\mu\nu}(q,r,p)$ the gauge coupling $g$, thus defining the
vertex ${\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,p)$, according to

$\displaystyle{\cal G}^{abc}_{\alpha\mu\nu}(q,r,p)$ $\displaystyle=$
$\displaystyle g{\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,p)\,,$
$\displaystyle{\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,p)$ $\displaystyle=$
$\displaystyle f^{abc}{\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,p)\,,$ (4)

where $f^{abc}$ denotes the structure constants of the group SU(3).

Note that both ${\cal G}^{abc}_{\alpha\mu\nu}(q,r,p)$ and
${\rm{I}}\\!\Gamma^{abc}_{\alpha\mu\nu}(q,r,p)$ will be employed throughout
this work; therefore, a clear diagrammatic distinction between them has been
introduced in Fig. 1. At tree level,

${\rm{I}}\\!\Gamma_{\\!0}^{\alpha\mu\nu}(q,r,p)=(q-r)^{\nu}g^{\alpha\mu}+(r-p)^{\alpha}g^{\mu\nu}+(p-q)^{\mu}g^{\nu\alpha}\,.$
(5)

Figure 1:  Diagrammatic conventions for the fully dressed three-gluon vertex
(left) and its tree-level counterpart (right).

The cornerstone of the SM is the observation that if the function
${\bf\Pi}(q)$ develops a pole at $q^{2}=0$ (“massless pole”), then, in that
kinematic limit, the propagator of the gauge boson saturates at a non-zero
value [12, 13], according to the sequence (Euclidean space)

$\lim_{q^{2}\to 0}{\bf\Pi}(q)=m^{2}/q^{2}\,\,\Longrightarrow\,\,\lim_{q^{2}\to
0}\,\Delta^{-1}(q)=\lim_{q^{2}\to
0}\,(q^{2}+m^{2})\,\,\Longrightarrow\,\,\Delta^{-1}(0)=m^{2}\,.$ (6)

This effect is interpreted as the dynamical generation of a mass, which is
identified with the positive residue of the pole.

In the absence of elementary scalar fields, massless poles may arise due to a
variety of reasons [98, 99, 100, 101]. In Yang-Mills theories, the pole
formation proceeds through the fusion of two gluons or of a ghost-antighost
pair into a color-carrying scalar bound state of vanishing mass [14, 15, 54,
11, 55, 56, 57, 58, 49, 59], to be denoted by $\Phi^{a}$. In the case of the
gluon fusion that we will consider throughout this work, the formation of this
bound state is controlled by a special BSE; details of the solutions obtained
are given in the following sections.

There are three main structures associated with this nonperturbative process,
which are diagrammatically depicted in Fig. 2:

(i) The effective vertex describing the interaction between $\Phi^{a}$ and two
gluons, denoted by

$B^{abc}_{\mu\nu}(q,r,p)=if^{abc}B_{\mu\nu}(q,r,p)\,.$ (7)

(ii) The propagator of the massless composite scalar, denoted by

$D^{ab}_{\Phi}(q)=\frac{i\delta^{ab}}{q^{2}}\,.$ (8)

(iii) The transition amplitude, $I^{ab}_{\alpha}(q)=\delta^{ab}I_{\alpha}(q)$,
connecting a gluon $A^{a}_{\alpha}$ with a scalar $\Phi^{b}$; Lorentz
invariance imposes that

$I_{\alpha}(q)=q_{\alpha}I(q)\,,$ (9)

where $I(q)$ is a scalar form factor.

Figure 2: Left: The effective vertex $B^{abc}_{\mu\nu}(q,r,p)$, with Lorentz,
color, and momentum conventions indicated. Center: The propagator of the
massless composite scalar, $D_{\Phi}^{ab}(q)$. Right: Gluon-scalar transition
amplitude, $I^{ab}_{\alpha}(q)$.

The emergence of items (i)-(iii) modifies profoundly the structure of the
four-gluon kernel, ${\cal
T}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})$, entering in the
skeleton expansion of the three-gluon vertex SDE, see Fig. 4. In particular,
as shown diagrammatically in Fig. 3, we have that

${\cal
T}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})={\mathcal{K}}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})+{\mathcal{M}}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})\,,$
(10)

where ${\mathcal{K}}$ denotes the regular, pole-free term, while
${\mathcal{M}}^{mnbc}$ is given by

${\mathcal{M}}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})=B^{mnx}_{\rho\sigma}(p_{1},p_{2},q)D^{xe}_{\Phi}(q)B^{ebc}_{\mu\nu}(p_{3},p_{4},-q)\,,$
(11)

with $q=p_{1}+p_{2}=-p_{3}-p_{4}$.

Figure 3: Decomposition of the four-gluon scattering kernel, ${\cal
T}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})$, into a regular part,
${\cal K}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})$, and a massless
pole part, ${\cal M}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})$,
according to Eq. (10).

We emphasize that, despite appearances, the kernel ${\cal T}$ is one-particle
irreducible, because the propagator $D_{\Phi}(q)$, present in the term
${\mathcal{M}}$, does not correspond to a fundamental field, but rather to a
composite excitation, representing a multitude of gluon interactions.

The way how the pole in ${\mathcal{M}}$ causes the saturation of the gluon
propagator is through the three-gluon vertex. In particular, consider the SDE
for the three-gluon vertex, whose diagrammatic representation is given in Fig.
4; note that only the diagrams relevant to our purposes are shown (first row).
The insertion of ${\mathcal{M}}^{mnbc}$ in this SDE generates a new term for
the three-gluon vertex (second row), to be denoted by
$V_{\alpha\mu\nu}(q,r,p)$. Specifically,
${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,p)$ assumes the general form111Due to
the Bose symmetry of ${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,p)$, poles appear
also in the channels carrying momenta $r$ or $p$, and are proportional to
$I_{\mu}(r)$ or $I_{\nu}(p)$. However, when these legs are internal to Feynman
diagrams, they get annihilated due to the contraction with the corresponding
Landau-gauge propagators.

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

where $\Gamma_{\alpha\mu\nu}(q,r,p)$ represents the pole-free component, while
$V_{\alpha\mu\nu}(q,r,p)$ is given by

$V_{\alpha\mu\nu}(q,r,p)=I_{\alpha}(q)\left(\frac{i}{q^{2}}\right)iB_{\mu\nu}(q,r,p)\,,$
(13)

as may be read off from Fig. 4. Then, due to Eq. (9), the term
$V_{\alpha\mu\nu}(q,r,p)$ becomes

$V_{\alpha\mu\nu}(q,r,p)=-\left(\frac{q_{\alpha}}{q^{2}}\right)I(q)B_{\mu\nu}(q,r,p)\,.$
(14)

Evidently, $V_{\alpha\mu\nu}(q,r,p)$ satisfies
$P^{\alpha^{\prime}\alpha}(q)V_{\alpha\mu\nu}(q,r,p)=0$; in that sense, the
poles are said to be “longitudinally coupled”.

Now, since the SDE controlling the momentum evolution of the gluon propagator
depends on the vertex ${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,p)$, the component
$V_{\alpha\mu\nu}(q,r,p)$ provides the massless pole required for the
activation of the SM. Specifically, focusing on the part of the gluon
propagator proportional to $q_{\mu}q_{\nu}$, we obtain the characteristic
diagram shown in Fig. 5, composed by the “square” of the transition amplitude
$I_{\alpha}(q)$. In particular, in the limit $q\to 0$, we have the central
result [55, 56, 57]

$m^{2}=g^{2}I^{2}\,,$ (15)

where we have defined the short-hand notation $I:=I(0)$.

Figure 4: First line: SDE for the three-gluon vertex. Second line: The pole
induced to the three-gluon vertex due to the ${\cal M}$ component of ${\cal
T}$ in Eq. (10) (see also Fig. 3). Figure 5:  Relation between the gluon mass
and the transition amplitude, $I_{\alpha}(q)$.

Let us finally emphasize that, in order to elucidate the novel aspects of the
method, in Figs. 4 and 5 we have omitted all diagrams containing ghosts as
internal lines or four-gluon vertices. This omission poses no problem for the
following two reasons. First, the omitted graphs do not interfere in any way
with the results derived in what follows; for a brief discussion, see Sec. X.
Second, we do not rely on the vertex SDE for an accurate description of the
form factor $\mathit{L}_{{sg}}$ [see Eq. (35)], appearing in some of the
central formulas; instead, we use the available lattice results for this
quantity [102, 103, 104, 105, 106, 107, 108, 109, 110].

##  III Setting the stage: Key results and main properties

In this section we present a collection of relations and properties that serve
as prerequisites for deriving the main results of this work. Note that, while
several of these items have already appeared in the literature, here we focus
on aspects that are especially pertinent for the task at hand.

###  III.1 An important limit

The kinematic limit relevant for the gluon mass generation is that of
vanishing momentum transfer, or, in terms of the propagator diagram in Fig. 5,
the limit $q\to 0$. It is therefore important to determine the behavior of the
effective vertex $B_{\mu\nu}(q,r,p)$ in that limit.

The general tensorial decomposition of this vertex is given by

$B_{\mu\nu}(q,r,p)=B_{1}\,g_{\mu\nu}+B_{2}\,r_{\mu}r_{\nu}+B_{3}\,p_{\mu}p_{\nu}+B_{4}\,r_{\mu}p_{\nu}+B_{5}\,p_{\mu}r_{\nu}\,,$
(16)

and $B_{i}:=B_{i}(q,r,p)$. The Bose symmetry of the full vertex
$B^{abc}_{\mu\nu}(q,r,p)$ under the exchange of the two incoming gluons
imposes that $B^{abc}_{\mu\nu}(q,r,p)=B^{acb}_{\nu\mu}(q,p,r)$, from which
follows that $B_{\mu\nu}(q,r,p)=-B_{\nu\mu}(q,p,r)$. Then, setting in this
last relation $q=0$ ($p=-r$), we obtain
$B_{\mu\nu}(0,r,-r)=-B_{\nu\mu}(0,r,-r)$. Given that

$B_{\mu\nu}(0,r,-r)=B_{1}(0,r,-r)\,g_{\mu\nu}+C_{1}(0,r,-r)r_{\mu}r_{\nu}\,,$
(17)

where $C_{1}:=B_{2}+B_{3}-B_{4}-B_{5}$, we conclude immediately that

$B_{1}(0,r,-r)=0=C_{1}(0,r,-r)\,\,\,\Longrightarrow\,\,B_{\mu\nu}(0,r,-r)=0\,.$
(18)

Let us then consider the Taylor expansion of $B_{\mu\nu}(q,r,p)$ around $q=0$,

$B_{\mu\nu}(0,r,-r)=q^{\alpha}{\cal B}_{\alpha\mu\nu}(0,r,-r)+\ldots\,,$ (19)

where we have introduced the shorthand notation

${\cal B}_{\alpha\mu\nu}(0,r,-r):=\left[\frac{\partial}{\partial
q^{\alpha}}B^{\mu\nu}(q,r,\,-r-q)\right]_{q=0}\,,$ (20)

accompanied by the diagrammatic representation shown in Fig. 6, and the
ellipsis denotes terms of higher order in $q$. When taking the limit indicated
in Eq. (20) we consider $r$ to be independent of $q$, in which case, due to
conservation of four-momentum, $p$ depends on $q$, since $p=-q-r$; therefore,
$\partial p^{\mu}/\partial q^{\alpha}=-g^{\mu}_{\alpha}$.

The next step is to contract ${\cal B}_{\alpha\mu\nu}(0,r,-r)$ by
$P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(r)$; this contraction appears
naturally in the diagrams where $B_{\mu\nu}(q,r,p)$ is inserted, because the
two gluon propagators attached to it are in the Landau gauge. Then, since the
tensorial decomposition of ${\cal B}_{\alpha\mu\nu}(0,r,-r)$ is given by

${\cal B}_{\alpha\mu\nu}(0,r,-r)={\cal B}_{1}(r)\,r_{\alpha}g_{\mu\nu}+{\cal
B}_{2}(r)\left[r_{\alpha}g_{\mu\nu}+r_{\mu}g_{\alpha\nu}+r_{\nu}g_{\mu\alpha}\right]+{\cal
B}_{3}(r)\,r_{\alpha}r_{\mu}r_{\nu}\,,$ (21)

only the first term survives this contraction, i.e.,

${\cal
B}_{\alpha\mu\nu}(0,r,-r)P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(r)={\cal
B}_{1}(r)\,r_{\alpha}P_{\mu^{\prime}\nu^{\prime}}(r)\,.$ (22)

It is now straightforward to establish that

$r_{\alpha}{\cal B}_{1}(r)=\left[\frac{\partial}{\partial
q^{\alpha}}B_{1}(q,r,\,-r-q)\right]_{q=0}=2r_{\alpha}\underbrace{\left[\frac{\partial{B}_{1}(q,r,p)}{\partial
p^{2}}\right]_{q=0}}_{{\mathbb{B}}(r)}\,,$ (23)

or

${\cal B}_{1}(r)=2{\mathbb{B}}(r)\,,$ (24)

from which follows that

${\cal
B}_{\alpha\mu\nu}(0,r,-r)=2\,\mathbb{B}(r)\,r_{\alpha}g_{\mu\nu}+\ldots\,,$
(25)

and then, from Eq. (19),

$B_{\mu\nu}(0,r,-r)=2\,(q\cdot r)\,\mathbb{B}(r)\,g_{\mu\nu}+\ldots\,,$ (26)

where the ellipses denote terms that get annihilated upon the aforementioned
contraction. The diagrammatic representation of Eq. (25) is shown in Fig. 6,
and will be used in the analysis that follows.

Figure 6: Diagrammatic representation of Eq. (25), corresponding to the first
nonvanishing term in the Taylor expansion of $B_{\mu\nu}^{abc}(q,r,p)$ around
$q=0$.

###  III.2 The gluon-scalar transition amplitude

Given the central relation of Eq. (15), we next derive a convenient expression
that allows us to compute the form factor $I$ in terms of the basic quantities
entering in its diagrammatic definition, see Fig. 2, item (iii).

For general momentum $q$, the transition amplitude ${I}^{\alpha}(q)$ is given
by

$I^{\alpha}(q)=-\frac{iC_{\rm
A}}{2}\int_{k}\Gamma_{\\!0}^{\alpha\beta\lambda}(q,k,-k-q)\Delta_{\beta\mu}(k)\Delta_{\lambda\nu}(k+q)B^{\mu\nu}(-q,-k,k+q)\,,$
(27)

where Eq. (7) was employed, the symmetry factor $\frac{1}{2}$ has been
supplied, and $C_{\rm A}$ is the Casimir eigenvalue of the adjoint
representation [$N$ for SU($N$)]. In addition, we use the shorthand notation

$\int_{k}:=\frac{1}{(2\pi)^{4}}\int\\!\\!{\rm d}^{4}k\,,$ (28)

where the use of a symmetry-preserving regularization scheme is implicitly
assumed.

Then, by virtue of Eq. (9), it is elementary to establish that

$I=\frac{1}{4}\left[\frac{\partial I^{\alpha}(q)}{\partial
q^{\alpha}}\right]_{q=0}\,.$ (29)

So, from Eq. (27) we obtain

$\displaystyle 4I=$ $\displaystyle\,-\frac{iC_{\rm
A}}{2}\int_{k}\left[\frac{\partial}{\partial
q^{\alpha}}\Gamma_{\\!0}^{\alpha\beta\lambda}(q,k,-k-q)\Delta_{\beta\mu}(k)\Delta_{\lambda\nu}(k+q)\right]_{q=0}\\!\\!\\!\\!\\!\\!B^{\mu\nu}(0,-k,k)$
$\displaystyle\,-\frac{iC_{\rm
A}}{2}\int_{k}\Gamma_{\\!0}^{\alpha\beta\lambda}(0,k,-k)\Delta_{\beta}^{\mu}(k)\Delta_{\lambda}^{\nu}(k)\,{\cal
B}_{\alpha\mu\nu}(0,-k,k)\,,$ (30)

and since $B_{\mu\nu}(0,-k,k)=0$ [see Eq. (18)], we have

$4I=-\frac{iC_{\rm
A}}{2}\int_{k}\Gamma_{\\!0}^{\alpha\beta\lambda}(0,k,-k)\Delta_{\beta}^{\mu}(k)\Delta_{\lambda}^{\nu}(k)\,{\cal
B}_{\alpha\mu\nu}(0,-k,k)\,.$ (31)

To further evaluate Eq. (30), use Eq. (5) to get

$\Gamma_{\\!0}^{\alpha\beta\lambda}(0,k,-k)=2k^{\alpha}g^{\beta\lambda}-k^{\lambda}g^{\alpha\beta}-k^{\beta}g^{\alpha\lambda}\,,$
(32)

and then Eq. (25) (with $r\to-k$), to finally obtain

$I=-\frac{3iC_{\rm A}}{2}\int_{k}k^{2}\Delta^{2}(k)\mathbb{B}(k)\,.$ (33)

Figure 7: Diagrammatic definition of the scalar form factor, $I$, of the
transition amplitude.

###  III.3 SDE of the three-gluon vertex and its soft gluon limit

The soft-gluon limit of the SDE for $\Gamma_{\alpha\mu\nu}^{abc}(q,r,p)$ is of
central importance for the ensuing analysis; indeed, as will be shown in
detail in Sec. VII, its proper use leads to considerable simplifications at
the level of the mass equation. In order to elucidate this important aspect,
we employ a simplified version of the SDE, which contains precisely the
structure required for certain crucial relations to be triggered.
Specifically, we consider the SDE represented in Fig. 8, composed only by the
tree-level contribution $(a_{1})$ and the diagram ($a_{2}$), while diagrams
where the gluon with momentum $q$ couples to a ghost-gluon or a four-gluon
vertex are omitted222The role of these terms will be discussed in section Sec.
X.. The reason why diagram ($a_{2}$) is singled out is because its kernel
${\cal K}$ appears also in the BSE for $\mathbb{B}(r)$, shown in Fig. 10.
This, in turn, will be instrumental for the implementation of the
renormalization procedure outlined in Sec. VII

Note that, in the standard version of this SDE, the three-gluon vertex in
diagram $(a_{2})$ is kept at tree level. The form employed here, with the
three-gluon vertex fully-dressed, corresponds to the BSE analogue of this
equation [111, 112, 113, 114, 115, 116]. Consequently, the skeleton expansion
of the kernel $\cal K$, given in Fig. 9, does not contain certain classes of
diagrams (e.g., ladder graphs) in order to avoid overcounting. In fact, ghost
loops aside, the diagrams of Fig. 9 comprise precisely the kernel of the
standard glueball BSE [117, 118]. The main advantage of this version of the
SDE is that the additional fully-dressed vertex absorbs the vertex
renormalization $Z_{3}$, defined in Eq. (44), which otherwise would be
multiplying the tree-level vertex; as a result, renormalization may be carried
out subtractively rather than multiplicatively. In particular, the validity of
the special formula in Eq. (61) hinges precisely on graph ($a_{2}$) having a
fully-dressed three-gluon vertex.

Figure 8: SDE for the regular part of the three-gluon vertex,
$\Gamma^{abc}_{\alpha\mu\nu}(q,r,p)$, in BS form.

We next factor out the color structure $f^{abc}$ by setting
$\Gamma_{\alpha\mu\nu}^{abc}(q,r,p)=f^{abc}\Gamma_{\alpha\mu\nu}(q,r,p)$, and
focus on the soft gluon limit, $\Gamma_{\alpha\mu\nu}(0,r,-r)$. The reason for
this choice is that, eventually, this vertex will be used in the graph of Eq.
(5), which is evaluated precisely in this kinematic limit. Setting $q=0$ into
Fig. 8, we get

$\Gamma_{\alpha\mu\nu}^{abc}(0,r,-r)=\Gamma_{\\!0\,\alpha\mu\nu}^{abc}(0,r,-r)\,+\,\int_{k}\Gamma^{axe}_{\alpha\gamma\delta}(0,k,-k)\Delta_{xm}^{\gamma\rho}(k)\Delta_{en}^{\delta\sigma}(k){\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)\,.$ (34)

Figure 9: Skeleton expansion of the scattering kernel, ${\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(p_{1},p_{2},p_{3},p_{4})$. The one-gluon exchange
is highlighted and the ellipsis denotes contributions with two or more loops.

To proceed further, consider the contraction of
$\Gamma_{\alpha\mu\nu}(0,r,-r)$ by
$P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(r)$; it is elementary to
establish that this vertex projection is described by a single form factor,
denoted by $\mathit{L}_{{sg}}(r)$, namely [106, 59]

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

Note that, in the Landau gauge that we employ, the
$\Gamma^{axe}_{\alpha\gamma\delta}(0,k,-k)$ on the r.h.s. of Eq. (34) is
automatically contracted by two transverse projectors, thus triggering Eq.
(35) therein. Then, contracting both sides of Eq. (34) by
$P^{\mu}_{\mu^{\prime}}(r)P^{\nu}_{\nu^{\prime}}(r)$, we obtain

$\displaystyle
f^{abc}\mathit{L}_{{sg}}(r)r_{\alpha}P_{\mu^{\prime}\nu^{\prime}}(r)=$
$\displaystyle\,f^{abc}r_{\alpha}P_{\mu^{\prime}\nu^{\prime}}(r)$ (36)
$\displaystyle-$
$\displaystyle\,f^{amn}\int_{k}k_{\alpha}\mathit{L}_{{sg}}(k)\Delta^{2}(k)P^{\rho\sigma}(k){\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)P^{\mu}_{\mu^{\prime}}(r)P^{\mu}_{\mu^{\prime}}(r)\,.$

Lastly, we contract Eq. (36) with
$f^{abc}r^{\alpha}g^{\mu^{\prime}\nu^{\prime}}$ to eliminate the color and
Lorentz indices. Using that $P^{\mu}_{\mu}(r)=3$, and that for SU($N$),
$f^{abc}f^{abc}=C_{\rm A}(N^{2}-1)$, we obtain

$\mathit{L}_{{sg}}(r)=\,1+\alpha_{s}\\!\int_{k}k^{2}\Delta^{2}(k)\,{K}(r,k)\,\mathit{L}_{{sg}}(k)\,,$
(37)

with the kernel $K(r,k)$ defined as

$\alpha_{s}K(r,k):=-\frac{(r\cdot
k)}{c\,r^{2}k^{2}}f^{abc}f^{amn}P^{\mu\nu}(r)P^{\rho\sigma}(k)\,{\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)\,,$ (38)

where we explicitly factor out the strong charge $\alpha_{s}=g^{2}/4\pi$. The
numerical factor $c:=3C_{\rm A}(N^{2}-1)$ arises from the projections of
Lorentz and color indices; for $N=3$, $c=72$.

###  III.4 Displacement of the Ward identity of the three-gluon vertex

In addition to triggering the SM, the massless poles are crucial for
preserving gauge invariance: their inclusion in the off-shell interaction
vertices guarantees that the Slavnov-Taylor identities (STIs) [119, 120] of
the theory retain the exact same form before and after mass generation [121].
In the case of the three-gluon vertex, the presence of these poles gives rise
to a characteristic effect, namely the displacement of the WI satisfied by the
soft-gluon form factor $\mathit{L}_{{sg}}(r)$ [59, 76], introduced in Eq.
(35).

In particular, let us denote by $\mathit{L}_{{sg}}(r)$ and
$\mathit{L}_{{sg}}^{\star}(r)$ the corresponding form factor when the SM is
turned on and off, respectively. When the SM is off, the three-gluon vertex is
simply given by $\Gamma_{\alpha\mu\nu}(q,r,p)$, and the STI it satisfies is
given by [4, 122, 123, 124]

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

where ${\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,p):=P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(p)$,
$F(q)$ denotes the ghost dressing function, defined from the ghost propagator,
$D^{ab}(q)$, by $D^{ab}(q)=i\delta^{ab}F(q)/q^{2}$, and $H_{\nu\mu}(p,q,r)$ is
the ghost-gluon kernel.

Taking the limit $q\to 0$ of both sides of Eq. (39), also known as the “soft-
gluon limit”, one obtains the corresponding WI, which simply states that

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

When the SM is activated, the vertex that satisfies the STI of Eq. (39) is the
full ${\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,k)$, given by Eq. (12). In
particular, with the aid of Eq. (14), the l.h.s. of Eq. (39) becomes

$q^{\alpha}{\rm{I}}\\!\Gamma_{\alpha\mu\nu}(q,r,p)\,{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,p)=[q^{\alpha}\Gamma_{\alpha\mu\nu}(q,r,p)-I(q)B_{\mu\nu}(q,r,p)]{\cal
T}_{\mu^{\prime}\nu^{\prime}}^{\mu\nu}(r,p)\,,$ (41)

while the r.h.s. remains exactly the same. Then, as has been shown in [59,
76], the limit $q\to 0$ of the STI yields

$\mathit{L}_{{sg}}(r^{2})=\mathit{L}_{{sg}}^{\\!\star}(r^{2})+{\mathbb{C}}(r^{2})\,,$
(42)

where

${\mathbb{C}}(r):=-I\,\mathbb{B}(r)\,,$ (43)

is denominated the “displacement function”.

The importance of this result is that one may determine ${\mathbb{C}}(r)$ as
the difference between the $\mathit{L}_{{sg}}(r^{2})$, extracted directly from
the lattice, and the expression for $\mathit{L}_{{sg}}^{\\!\star}(r^{2})$ on
the r.h.s. of Eq. (40), which is evaluated using lattice inputs for its
various components [59, 76]. The result for ${\mathbb{C}}(r)$ that emerges
through this procedure will serve as benchmark for the numerical analysis
presented in Sec. IX.

##  IV Renormalization: general considerations

In this section we discuss in detail the renormalization of the key dynamical
components that emerge from the activation of the SM. In particular, we will
focus on the renormalization properties of the central quantities $I$ and
$B_{\mu\nu}$.

We focus exclusively on the gluon propagator and the three-gluon vertex, which
are the two basic ingredients comprising the structures considered in this
work. Denoting by the index “R” the renormalized quantities, we have

$\Delta_{{\scriptscriptstyle
R}}(q^{2})=Z^{-1}_{A}\Delta(q^{2})\,,\quad\quad\quad{\rm{I}}\\!\Gamma^{\alpha\mu\nu}_{\\!\\!{\scriptscriptstyle
R}}(q,r,p)=Z_{3}{\rm{I}}\\!\Gamma^{\alpha\mu\nu}(q,r,p)\,,\quad\quad\quad
g_{{\scriptscriptstyle R}}=Z_{g}^{-1}g\,,$ (44)

where $Z^{1/2}_{A}$ is the wave function renormalization constants of the
gluon field, $Z_{3}$ the renormalization constant of the three-gluon vertex,
and $Z_{g}$ the coupling renormalization constant, given by

$Z_{g}=Z_{3}Z_{A}^{-3/2}\,.$ (45)

Furthermore, we will introduce two additional renormalization constants, to be
denoted by $Z_{I}$ and $Z_{B}$, which renormalize the quantities $I$ and
$B^{\mu\nu}(q,r,p)$, respectively, according to [125]

$I_{{\scriptscriptstyle R}}=Z^{-1}_{I}I\,,\quad\quad B_{{\scriptscriptstyle
R}}^{\mu\nu}(q,r,p)=Z^{-1}_{B}B^{\mu\nu}(q,r,p)\,.$ (46)

The partial derivatives of $B^{\mu\nu}(q,r,p)$ also renormalize in the same
way, namely

${\cal B}_{{\scriptscriptstyle R}}^{\alpha\mu\nu}(0,r,-r)=Z^{-1}_{B}{\cal
B}_{{\scriptscriptstyle
R}}^{\alpha\mu\nu}(0,r,-r)\,,\quad\quad\mathbb{B}_{{\scriptscriptstyle
R}}(r)=Z^{-1}_{B}\mathbb{B}(r)\,.$ (47)

Given that both the $I$ and $B^{\mu\nu}$ are composed by the fundamental
fields and vertices, it is natural to expect that the $Z_{I}$ and $Z_{B}$
should be expressed in terms of the $Z_{A}$ and $Z_{3}$ introduced in Eq.
(44); indeed, we find that [125]

$Z_{I}=Z^{-1}_{3}Z_{A}\,,\quad\quad\quad Z_{B}=Z^{-1}_{A}\,.$ (48)

To see how the relations in Eq. (48) arise, we first turn to Eq. (15). Since
$m^{2}:=\Delta^{-1}(0)$, from the first relation in Eq. (44) we have that
[125, 57]

$m^{2}=Z^{-1}_{A}m^{2}_{{\scriptscriptstyle R}}\,.$ (49)

On the other hand, combining Eq. (15), the first relation in Eq. (46), and Eq.
(45), we have

$m^{2}=g^{2}I^{2}=Z^{2}_{g}\,Z^{2}_{I}\,\underbrace{g^{2}_{{\scriptscriptstyle
R}}I^{2}_{{\scriptscriptstyle R}}}_{m^{2}_{{\scriptscriptstyle
R}}}=Z^{2}_{3}\,Z^{-3}_{A}\,Z^{2}_{I}\,m^{2}_{{\scriptscriptstyle R}}\,.$ (50)

Then, the direct comparison of Eqs. (49) and (50) leads immediately to the
first relation of Eq. (48).

The second relation in Eq. (48) may be obtained from Eq. (43), by noticing
that, since $V_{\alpha\mu\nu}(q,r,p)$ is a component of the three-gluon vertex
${\rm{I}}\\!\Gamma^{\alpha\mu\nu}(q,r,p)$, it is renormalized as in Eq. (44),
namely

$V^{\alpha\mu\nu}(q,r,p)=Z^{-1}_{3}\,V^{\alpha\mu\nu}_{\\!\\!{\scriptscriptstyle
R}}(q,r,p)\,,$ (51)

and, therefore,

${\mathbb{C}}(k)=Z^{-1}_{3}\,{\mathbb{C}}_{\\!{\scriptscriptstyle R}}(k)\,.$
(52)

On the other hand, from Eq. (43), using Eq. (46) and the first relation of Eq.
(48), we have that

${\mathbb{C}}(k)=-Z_{I}Z_{B}\underbrace{I_{{\scriptscriptstyle
R}}\mathbb{B}_{{\scriptscriptstyle
R}}(k)}_{-{\mathbb{C}}_{\\!{\scriptscriptstyle
R}}(k)}=Z^{-1}_{3}Z_{A}Z_{B}\,{\mathbb{C}}_{\\!{\scriptscriptstyle R}}(k)\,.$
(53)

Then, the comparison between Eqs. (52) and (53) yields the second relation of
Eq. (48).

We next turn to the renormalization of the kernel ${\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(q,r,p,t)$, appearing in Eq. (34), and later on in
Eq. (65). To that end, notice that ${\cal K}$ is a part of the four-gluon
amplitude ${\cal G}^{mnbc}_{\rho\sigma\mu\nu}(q,r,p,t)=\langle
0|T[\widetilde{A}^{m}_{\rho}(q)\widetilde{A}^{n}_{\sigma}(r)\widetilde{A}^{b}_{\mu}(p)\widetilde{A}^{c}_{\nu}(t)]|0\rangle$,
with the external legs amputated, i.e.,

${\cal G}(q,r,p,t)=\Delta(q)\Delta(r)\Delta(p)\Delta(t){\cal
K}(q,r,p,t)+\cdots$ (54)

where the ellipsis denotes the diagrams excluded when switching from the SDE
to the BSE kernel, as discussed in Sec. III.3. Since,
$A^{a\mu}_{{\scriptscriptstyle R}}=Z^{1/2}_{A}A^{a\mu}$, the definition of
${\cal G}$ in terms of gauge fields implies that ${\cal
G}_{{\scriptscriptstyle R}}=Z^{-2}_{A}\,{\cal G}$. Therefore, from Eqs. (54)
and (44), we get

${\cal K}_{{\scriptscriptstyle R}}(q,r,p,t)=Z^{2}_{A}\,{\cal K}(q,r,p,t)\,,$
(55)

and, since ${\cal K}\sim\alpha_{s}K$, we have

${K}_{{\scriptscriptstyle R}}(q,r,p,t)=Z^{2}_{A}Z^{2}_{g}\,{K}(q,r,p,t)\,.$
(56)

From the above relations it is immediate to establish that the combinations
$\Delta{\mathbb{B}}$, $\Delta^{2}{\cal K}$, and $\alpha_{s}\Delta^{2}K$ are
renormalization-group invariant (RGI), i.e.,

$\Delta\mathbb{B}=\Delta_{{\scriptscriptstyle
R}}\mathbb{B}_{{\scriptscriptstyle R}}\,,\qquad\qquad\Delta^{2}{\cal
K}=\Delta^{2}_{{\scriptscriptstyle R}}\,{\cal K}_{{\scriptscriptstyle
R}}\,,\qquad\qquad\alpha_{s}\Delta^{2}K=\alpha_{s}^{{\scriptscriptstyle
R}}\Delta^{2}_{{\scriptscriptstyle R}}K_{{\scriptscriptstyle R}}\,.$ (57)

Note that, by virtue of the renormalization rule $B_{{\scriptscriptstyle
R}}^{\mu\nu}=Z_{A}B^{\mu\nu}$, given by Eqs. (46) and (48), the kernel
${\mathcal{M}}$ defined in Eq. (11) renormalizes as ${\cal
M}_{{\scriptscriptstyle R}}=Z^{2}_{A}\,{\cal M}$, i.e., , exactly as the
kernel ${\mathcal{K}}$ in Eq. (55); this is consistent with the fact that both
${\mathcal{K}}$ and ${\mathcal{M}}$ are parts of the same four-gluon kernel,
see Eq. (10) and Fig. 4.

Armed with the above results, we may now derive a useful expression for the
$I_{{\scriptscriptstyle R}}$, using Eq. (31) as our point of departure.
Substituting the bare quantities comprising Eq. (31) by renormalized ones, we
have

$4Z_{I}I_{{\scriptscriptstyle R}}=-\frac{iC_{\rm
A}}{2}Z^{2}_{A}Z_{B}\int_{k}\Gamma_{\\!0}^{\alpha\beta\lambda}(0,k,-k)\Delta_{{{\scriptscriptstyle
R}}\,\beta}^{\mu}(k)\Delta_{{{\scriptscriptstyle R}}\,\lambda}^{\nu}(k)\,{\cal
B}_{\\!{{\scriptscriptstyle R}}\,\alpha\mu\nu}(0,-k,k)\,,$ (58)

and, after employing Eq. (48),

$4I_{{\scriptscriptstyle R}}=-\frac{iC_{\rm
A}}{2}Z_{3}\int_{k}\Gamma_{\\!0}^{\alpha\beta\lambda}(0,k,-k)\Delta_{{{\scriptscriptstyle
R}}\,\beta}^{\mu}(k)\Delta_{{{\scriptscriptstyle R}}\,\lambda}^{\nu}(k)\,{\cal
B}_{\\!{{\scriptscriptstyle R}}\,\alpha\mu\nu}(0,-k,k)\,,$ (59)

or, from Eq. (33),

$I_{{\scriptscriptstyle R}}=-\frac{3iC_{\rm
A}}{2}Z_{3}\int_{k}k^{2}\Delta_{{\scriptscriptstyle
R}}^{2}(k)\mathbb{B}_{{\scriptscriptstyle R}}(k)\,.$ (60)

We finally turn to the renormalization of Eq. (37). By virtue of Eq. (35), it
is clear that $\mathit{L}_{{sg}}^{\\!{\scriptscriptstyle
R}}(r)=Z_{3}\mathit{L}_{{sg}}(r)$. Then, using Eq. (57), it is straightforward
to show that the renormalized version of Eq. (37) is given by

$\mathit{L}_{{sg}}^{\\!{\scriptscriptstyle
R}}(r)=\,Z_{3}+\alpha_{s}^{{\scriptscriptstyle
R}}\\!\\!\int_{k}k^{2}\Delta^{2}_{{\scriptscriptstyle
R}}(k)\,{K}_{{\scriptscriptstyle
R}}(r,k)\,\mathit{L}_{{sg}}^{\\!{\scriptscriptstyle R}}(k)\,.$ (61)

As announced in Sec. III.3, the renormalization required for Eq. (61) is
subtractive.

In what follows we will suppress the index “R” in order to avoid notational
clutter.

##  V Non-linear Bethe-Salpeter equation for pole formation

In this section we derive the BSE satisfied by the amplitude $\mathbb{B}(r)$.
This type of equation has been derived in the literature before [55, 56, 57,
58, 59]; however, in contradistinction to these earlier versions, the present
BSE is not linearized, having its cubic nature fully retained.

The diagrammatic form of the BSE for general $q$ is shown in the left part of
Fig. 10, where the four-gluon kernel ${\cal T}$ is depicted in Fig. 3.
Specifically, from the first equality of Fig. 10 we have

$B^{abc}_{\mu\nu}(q,r,p)=(G_{\cal T})^{abc}_{\mu\nu}(q,r,p)\,,$ (62)

where

$(G_{\cal
T})^{abc}_{\mu\nu}(q,r,p)=\int_{k}B^{axe}_{\alpha\beta}(q,k,-k-q)\Delta_{xm}^{\alpha\rho}(k)\Delta_{en}^{\beta\sigma}(k+q){\cal
T}^{mnbc}_{\rho\sigma\mu\nu}(-k,k+q,r,p)\,.$ (63)

Then, after using Eqs. (10) and (11), two distinct terms arise, namely

$(G_{\cal T})^{abc}_{\mu\nu}(q,r,p)=(G_{\cal
K})^{abc}_{\mu\nu}(q,r,p)+(G_{\\!{\cal M}})^{abc}_{\mu\nu}(q,r,p)\,,$ (64)

with

$\displaystyle(G_{\cal K})^{abc}_{\mu\nu}(q,r,p)$ $\displaystyle=$
$\displaystyle\int_{k}B^{axe}_{\alpha\beta}(q,k,-k-q)\Delta_{xm}^{\alpha\rho}(k)\Delta_{en}^{\beta\sigma}(k+q){\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k+q,r,p)\,,$ $\displaystyle(G_{\\!{\cal
M}})^{abc}_{\mu\nu}(q,r,p)$ $\displaystyle=$
$\displaystyle\Omega^{ad}(q)D^{ds}_{\Phi}(q)B^{sbc}_{\mu\nu}(q,r,p)\,,$ (65)

where (symmetry factor $\frac{1}{2}$ included)

$\Omega^{ad}(q)=\frac{1}{2}\int_{k}B^{axe}_{\alpha\beta}(q,k,-k-q)\Delta_{xm}^{\alpha\rho}(k)\Delta_{en}^{\beta\sigma}(k+q)B^{dnm}_{\sigma\rho}(-q,k+q,-k)\,.$
(66)

Note that, while the term $(G_{\cal K})$ is linear in $B^{abc}_{\mu\nu}$, the
term $(G_{\\!{\cal M}})$ is cubic. As a result, the scale ambiguity of the
solutions, present when only the term $(G_{\cal K})$ is considered, is now
eliminated (see next section).

Next, it is important to establish that the product
$\Omega^{ad}(q)D^{ds}_{\Phi}(q)$ appearing in $(G_{\\!{\cal M}})$ is finite as
$q\to 0$; this is so, because, in that limit, $\Omega^{ad}(q)\sim
q^{2}\delta^{ad}$, thus cancelling exactly the massless pole in
$D^{ds}_{\Phi}(q)$.

To derive this key result from Eq. (66), set first
$\Omega^{ad}(q)=\delta^{ad}\Omega(q)$ and carry out the color algebra to
obtain

$\Omega(q)=\frac{C_{\rm
A}}{2}\int_{k}B_{\alpha\beta}(q,k,-k-q)P^{\alpha\rho}(k)P^{\beta\sigma}(k+q)B_{\rho\sigma}(-q,-k,k+q)\Delta(k)\Delta(k+q)\,.$
(67)

Then, taking the limit $q\to 0$ of Eq. (67) using Eq. (26), we find

$\displaystyle\lim_{q\to 0}\Omega(q)$ $\displaystyle=$
$\displaystyle\frac{C_{\rm
A}}{2}\int_{k}{B}_{\alpha\beta}(0,k,-k)P^{\alpha\rho}(k)P^{\beta\sigma}(k){B}_{\sigma\rho}(0,k,-k)\Delta^{2}(k)$
(68) $\displaystyle=$ $\displaystyle 6C_{\rm A}\int_{k}(q\cdot
k)^{2}\Delta^{2}(k){\mathbb{B}}^{2}(k)\,,$

and thus

$\lim_{q\to
0}\Omega(q)=q^{2}{\widetilde{\omega}}\,,\qquad\qquad{\widetilde{\omega}}:=\frac{3C_{\rm
A}}{2}\int_{k}k^{2}\Delta^{2}(k){\mathbb{B}}^{2}(k)\,.$ (69)

With this result in hand, it is immediate to establish that

$(G_{\\!{\cal
M}})^{abc}_{\mu\nu}(0,r,-r)=\omega\,{B}_{\mu\nu}^{abc}(0,r,-r)\,,\qquad\omega:=i\,{\widetilde{\omega}}\,;$
(70)

therefore, this term combines directly with the term $B^{abc}_{\mu\nu}(q,r,p)$
on the l.h.s of Eq. (62) in the same limit. In particular, after defining the
new variable,

$t:=1-\omega\,,$ (71)

Eq. (62) becomes

$t\,{B}^{abc}_{\mu\nu}(0,r,-r)=\int_{k}B^{axe}_{\alpha\beta}(0,k,-k)\Delta_{xm}^{\alpha\rho}(k)\Delta_{en}^{\beta\sigma}(k){\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)\,.$ (72)

Given Eq. (18), it is clear that, in the limit $q=0$, Eq. (72) yields to
lowest order a trivial result ($0=0$). So, in order to extract nontrivial
information from it, one must equate the terms linear in $q$ on both of its
sides. Specifically, using Eq. (19) for the $B_{\mu\nu}(0,r-r)$ and
$B_{\alpha\beta}(0,k-k)$ appearing in Eq. (72), one obtains the equation333We
assume that $t\neq 0$; if $t=0$, the l.h.s. of Eq. (72) vanishes, and the
resulting equation is no longer of the BSE type.

${\cal B}^{abc}_{\lambda\mu\nu}(0,r,-r)=t^{-1}\int_{k}{\cal
B}^{axe}_{\lambda\alpha\beta}(0,k,-k)\Delta_{xm}^{\alpha\rho}(k)\Delta_{en}^{\beta\sigma}(k){\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)\,,$ (73)

which admits the diagrammatic representation given in the second line of Fig.
10.

Then, contracting both sides by
$P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(r)$ and employing Eq. (22),
we get

$f^{abc}\,\mathbb{B}(r)\,r_{\lambda}P_{\mu^{\prime}\nu^{\prime}}(r)\,=-t^{-1}f^{amn}\int_{k}\mathbb{B}(k)\,\Delta^{2}(k)\,k_{\lambda}P^{\sigma\rho}(k)\,{\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)P_{\mu^{\prime}}^{\mu}(r)P_{\nu^{\prime}}^{\nu}(r)\,,$
(74)

and, after contracting both sides by
$f^{abc}\,r^{\lambda}g^{\mu^{\prime}\nu^{\prime}}$, and using that
$P_{\mu}^{\mu}(r)=3$, we arrive at the following equation for the amplitude
$\mathbb{B}(r)$,

$\mathbb{B}(r)=t^{-1}\alpha_{s}\int_{k}k^{2}\Delta^{2}(k){K}(r,k)\mathbb{B}(k)\,,$
(75)

where $K(r,k)$ is exactly the kernel defined in Eq. (38).

We end this section by pointing out that Eq. (75) is RGI. This follows
directly by observing that, due to the first relation in Eq. (57),
${\widetilde{\omega}}$ itself, defined through Eq. (69), is RGI, and so is the
$t$ in Eq. (71) i.e., $\omega=\omega_{{\scriptscriptstyle R}}$ and
$t=t_{{\scriptscriptstyle R}}$. Then, the full BSE in Eq. (75) is RGI due to
the third relation in Eq. (57), and the fact that the $Z_{B}$ introduced to
renormalize the $\mathbb{B}$ cancels on both sides. As a result, one may
substitute into Eq. (75) directly bare for renormalized quantities, without
any renormalization constants appearing in the final expression, i.e.,

$\mathbb{B}_{{\scriptscriptstyle R}}(r)=t_{{\scriptscriptstyle
R}}^{-1}\alpha_{s}^{{\scriptscriptstyle
R}}\int_{k}k^{2}\Delta_{{\scriptscriptstyle R}}^{2}(k){K}_{{\scriptscriptstyle
R}}(r,k)\mathbb{B}_{{\scriptscriptstyle R}}(k)\,.$ (76)

Again, the index “R” will be suppressed in what follows.

Figure 10: BSE for the effective vertex $B^{abc}_{\mu\nu}(q,r,p)$ (left) and
its $q=0$ limit (right).

##  VI Dynamical Scale-fixing

In order to proceed further with our analysis, it is necessary to pass certain
key equations to Euclidean space, following standard transformation rules.
Specifically, we consider all physical momenta to be space-like; in
particular, $r^{2}\to-r_{{\rm{\scriptscriptstyle E}}}^{2}$, where
$r_{{\rm{\scriptscriptstyle E}}}^{2}>0$. In addition, we use the conversion
relations

$\displaystyle\int_{k}=$ $\displaystyle\,i\int_{k_{{\rm{\scriptscriptstyle
E}}}}\,,\qquad$ $\displaystyle\Delta_{{\rm{\scriptscriptstyle
E}}}(r_{{\rm{\scriptscriptstyle E}}}^{2})=$
$\displaystyle\,-\Delta(-r_{{\rm{\scriptscriptstyle E}}}^{2})\,,\qquad$
$\displaystyle\mathbb{B}_{{\rm{\scriptscriptstyle
E}}}(r_{{\rm{\scriptscriptstyle E}}}^{2})=$
$\displaystyle\,-\mathbb{B}(-r^{2}_{{\rm{\scriptscriptstyle E}}})\,,$
$\displaystyle{\mathbb{C}}_{{\rm{\scriptscriptstyle
E}}}(r_{{\rm{\scriptscriptstyle E}}}^{2})=$
$\displaystyle\,-{\mathbb{C}}(-r^{2}_{{\rm{\scriptscriptstyle E}}})\,,\qquad$
$\displaystyle\mathit{L}_{{sg}}^{\\!{{\rm{\scriptscriptstyle
E}}}}(r_{{\rm{\scriptscriptstyle E}}}^{2})=$
$\displaystyle\,\mathit{L}_{{sg}}(-r_{{\rm{\scriptscriptstyle E}}}^{2})\,.$
(77)

Applying the above relations to Eqs. (60), (69), and (70), the expressions for
$\omega_{{\rm{\scriptscriptstyle E}}}$ and $I_{{\rm{\scriptscriptstyle E}}}$
read

$\displaystyle\omega_{{\rm{\scriptscriptstyle E}}}$ $\displaystyle=$
$\displaystyle\,\frac{3C_{\rm A}}{2}\int_{k_{{\rm{\scriptscriptstyle
E}}}}k_{{\rm{\scriptscriptstyle E}}}^{2}\Delta^{2}_{{\rm{\scriptscriptstyle
E}}}(k_{{\rm{\scriptscriptstyle E}}}^{2})\mathbb{B}_{{\rm{\scriptscriptstyle
E}}}^{2}(k_{{\rm{\scriptscriptstyle E}}}^{2})\,,$ (78) $\displaystyle
I_{{\rm{\scriptscriptstyle E}}}$ $\displaystyle=$
$\displaystyle\,\frac{3C_{\rm A}Z_{3}}{2}\int_{k_{{\rm{\scriptscriptstyle
E}}}}k_{{\rm{\scriptscriptstyle E}}}^{2}\Delta^{2}_{{\rm{\scriptscriptstyle
E}}}(k_{{\rm{\scriptscriptstyle E}}}^{2})\mathbb{B}_{{\rm{\scriptscriptstyle
E}}}(k_{{\rm{\scriptscriptstyle E}}}^{2})\,.$ (79)

The Euclidean versions of Eqs. (61) and (76) may be derived in a similar way.
In doing so, note that the kernel $K$ must be cast in the form444The origin of
the factor “$i$” may be easily understood at the level of the one-gluon
exchange diagram of Fig. 9; it is included in the definition of the gluon
propagator, see Eq. (1).

$K(r,k)=iK_{{\rm{\scriptscriptstyle E}}}(r_{{\rm{\scriptscriptstyle
E}}},k_{{\rm{\scriptscriptstyle E}}})\,,$ (80)

where

$i\,\alpha_{s}K_{{\rm{\scriptscriptstyle E}}}(r_{{\rm{\scriptscriptstyle
E}}},k_{{\rm{\scriptscriptstyle E}}}):=\frac{(r_{{\rm{\scriptscriptstyle
E}}}\cdot k_{{\rm{\scriptscriptstyle E}}})}{c\,r^{2}_{{\rm{\scriptscriptstyle
E}}}k^{2}_{{\rm{\scriptscriptstyle
E}}}}f^{abc}f^{amn}\left[P^{\mu\nu}(r)P^{\rho\sigma}(k)\,{\cal
K}^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)\right]_{{\rm{\scriptscriptstyle
E}}}\,,$ (81)

and $\left[\ldots\right]_{{\rm{\scriptscriptstyle E}}}$ denotes the
application of Eq. (77) to the expressions in square brackets.

Then, with $t_{{\rm{\scriptscriptstyle E}}}:=1-\omega_{{\rm{\scriptscriptstyle
E}}}$,

$\mathbb{B}_{{\rm{\scriptscriptstyle E}}}(r_{{\rm{\scriptscriptstyle
E}}}^{2})=t_{{\rm{\scriptscriptstyle
E}}}^{-1}\alpha_{s}\int_{k_{{\rm{\scriptscriptstyle
E}}}}k_{{\rm{\scriptscriptstyle E}}}^{2}\Delta_{{\rm{\scriptscriptstyle
E}}}^{2}(k_{{\rm{\scriptscriptstyle E}}}^{2}){K}_{{\rm{\scriptscriptstyle
E}}}(r_{{\rm{\scriptscriptstyle E}}},k_{{\rm{\scriptscriptstyle
E}}})\mathbb{B}_{{\rm{\scriptscriptstyle E}}}(k_{{\rm{\scriptscriptstyle
E}}}^{2})\,,$ (82)

whereas the equation for $\mathit{L}_{{sg}}^{{\rm{\scriptscriptstyle E}}}(r)$
reads

$\mathit{L}_{{sg}}^{\\!{{\rm{\scriptscriptstyle
E}}}}(r_{{\rm{\scriptscriptstyle
E}}}^{2})=\,Z_{3}+\alpha_{s}\\!\int_{k_{{\rm{\scriptscriptstyle
E}}}}k_{{\rm{\scriptscriptstyle E}}}^{2}\Delta_{{\rm{\scriptscriptstyle
E}}}^{2}(k_{{\rm{\scriptscriptstyle E}}}^{2})\,{K}_{{\rm{\scriptscriptstyle
E}}}(r_{{\rm{\scriptscriptstyle E}}},k_{{\rm{\scriptscriptstyle
E}}})\,\mathit{L}_{{sg}}^{\\!{{\rm{\scriptscriptstyle
E}}}}(k_{{\rm{\scriptscriptstyle E}}}^{2})\,.$ (83)

Finally, the expressions for the Euclidean displacement function and mass,
retain the same form as in Eqs. (15) and (43), but now carry indices “E”,
i.e.,

${\mathbb{C}}_{{\rm{\scriptscriptstyle E}}}(r_{{\rm{\scriptscriptstyle
E}}}^{2}):=-I_{{\rm{\scriptscriptstyle
E}}}\,\mathbb{B}_{{\rm{\scriptscriptstyle E}}}(r_{{\rm{\scriptscriptstyle
E}}}^{2})\,,\qquad\quad m^{2}=g^{2}I^{2}_{{\rm{\scriptscriptstyle E}}}\,.$
(84)

In order to simplify the notation, in what follows the index “E” will be
omitted.

We next focus on the BSE of Eq. (82). Strictly speaking, Eq. (82) is a non-
linear integral equation, due to the quadratic dependence of $t$ on the
unknown function $\mathbb{B}$. However, the fact that $t$ is a constant,
independent of the variable $r$, allows one to solve Eq. (82) as an eigenvalue
problem, typical of a linear homogeneous equation, with one crucial
difference: the solutions do not suffer from the scale indeterminacy intrinsic
to such a linear equation, because the presence of $t$ fixes the scale, up to
an overall sign.

To appreciate the role played by the parameter $t$, we temporarily set $t=1$.
Then, let us suppose that the kernel $K$ is such that the eigenvalue that
yields a nontrivial solution for $\mathbb{B}(k)$ requires that
$\alpha_{s}\to\alpha_{\star}$, where $\alpha_{\star}$ differs from the value
predicted for $\alpha_{s}$ within the renormalization scheme employed.

Then, restoring $t$ at the level of Eq. (76), it becomes clear that it has to
compensate precisely for the difference between $\alpha_{s}$ and
$\alpha_{\star}$. In particular, we must have

$t=\frac{\alpha_{s}}{\alpha_{\star}}\,,$ (85)

which brings Eq. (76) into the required form

$\mathbb{B}(r)=\alpha_{\star}\int_{k}k^{2}\Delta^{2}(k){K}(r,k)\mathbb{B}(k)\,.$
(86)

To understand how the presence of the parameter $t$ in Eq. (76) leads to the
determination of the scale of $\mathbb{B}(k)$, note that, by combining Eqs.
(71) and (85) we obtain the condition

${\omega}=1-\alpha_{s}/\alpha_{\star}\,.$ (87)

Since the value of the l.h.s. of Eq. (87) is fixed, ${\omega}$ is completely
determined, and, therefore, the size of the $\mathbb{B}(k)$ entering in the
definition of Eq. (78) is constrained.

To see this, let $\mathbb{B}_{0}(k)$ be a solution of Eq. (76) with a scale
set arbitrarily, e.g., by imposing that the global maximum of
$\mathbb{B}_{0}(k)$ is at $1$, and denote by $\omega_{0}$ the value of
$\omega$ when $\mathbb{B}(k)\to\mathbb{B}_{0}(k)$ is substituted in Eq. (78).
The $\mathbb{B}$ that is compatible with Eq. (87) is related to
$\mathbb{B}_{0}$ by a multiplicative constant, $\sigma$, i.e.,
$\mathbb{B}(k)=\sigma\,\mathbb{B}_{0}(k)$, whose value is determined from the
equation

$\sigma^{2}\omega_{0}=1-\alpha_{s}/\alpha_{\star}\,,$ (88)

or, equivalently,

$\sigma=\pm\sqrt{\frac{1-\alpha_{s}/\alpha_{\star}}{{{\omega}}_{0}}}\,.$ (89)

Note that, since ${{\omega}}$ is quadratic in $\mathbb{B}(k)$, the sign of
$\sigma$, and hence the sign of $\mathbb{B}$ itself, is left undetermined.
However, as we can see from Eqs. (79) and (84), ${\mathbb{C}}(k)$ is quadratic
in $\mathbb{B}(k)$, and therefore does not get affected by the sign ambiguity
of Eq. (89). In fact, the overall sign of ${\mathbb{C}}(k)$ turns out to be
negative for the entire range of Euclidean momenta, in agreement with the sign
found in the lattice extraction of ${\mathbb{C}}(k)$ presented in [59, 76].

##  VII Renormalization of the gluon mass equation

The compact relation in Eq. (15) expresses the gluon mass in terms of the
gluon-scalar transition amplitude $I$, which depends, in turn, on the BS
amplitude $\mathbb{B}(r)$. After the renormalization and the rotation to the
Euclidean space, $I$ is given by Eq. (79). The obvious difficulty associated
with the use of the latter equation is the renormalization constant $Z_{3}$
multiplying the relevant integral.

It turns out that a subtle sequence of relations allows one to implement the
multiplicative renormalization implicit in Eq. (79) without having to carry
out any integrations. The upshot of these considerations is that the
renormalization of $I$ finally amounts to the effective replacement
$Z_{3}\to\omega\mathit{L}_{{sg}}(k)$ at the level of Eq. (79), i.e.,

$I=\frac{3C_{\rm
A}}{2}\omega\int_{k}k^{2}\Delta^{2}(k)\mathit{L}_{{sg}}(k)\mathbb{B}(k)\,.$
(90)

In what follows we will present a detailed derivation of the above key result.

The first observation is that, by virtue of Eq. (61), the $Z_{3}$ in Eq. (60)
may be substituted by the combination

$Z_{3}=\,\mathit{L}_{{sg}}(k)-\alpha_{s}\\!\\!\int_{\ell}\ell^{2}\Delta^{2}(\ell)\,{K}(k,\ell)\,\mathit{L}_{{sg}}(\ell)\,.$
(91)

This procedure is typical when dealing with multiplicative renormalizability
at the level of the SDEs, see, e.g., [126], and is intimately connected with
the so-called “skeleton expansion” of the SDE kernel.

The second observation is novel, and completely specific to the particular
form of the BSE in Eq. (82). In particular, when the substitution of Eq. (91)
is implemented, the BSE of Eq. (82) is triggered inside Eq. (79), leading to a
crucial cancellation, and, finally, to Eq. (90).

In order to convey the underlying idea with a simple example, consider two
functions $g(x)$ and $f(x)$, satisfying the system of integral equations

$g(x)=z+\int\\!\\!dy\,\tilde{g}(y)K(x,y)\,,\qquad\qquad
f(x)=b^{-1}\\!\int\\!\\!dy\,\tilde{f}(y)K(x,y)\,,$ (92)

where $\tilde{g}(y):=g(y)\,u(y)$ and $\tilde{f}(y):=f(y)\,u(y)$, with $u(y)$ a
well-behaved function. The parameters $z$ and $b$ are real numbers, with
$b\neq 0,1$, and the limits of integration are arbitrary. Finally, the common
kernel $K(x,y)$ satisfies the symmetry relation $K(x,y)=K(y,x)$.

Let us further assume that the value of a constant $\beta$ is given by the
integral

$\beta=z\\!\int\\!\\!dx\,\tilde{f}(x)\,.$ (93)

The dependence of $\beta$ on the parameter $z$ may be eliminated in favor of
the function $g(x)$ by appealing to the system of Eq. (92). In particular, one
employs the following sequence of steps

$\displaystyle\beta$ $\displaystyle=$
$\displaystyle\int\\!\\!dx\,z\,\tilde{f}(x)$ (94) $\displaystyle=$
$\displaystyle\int\\!\\!dx\left[g(x)-\int\\!\\!dy\,\tilde{g}(y)K(x,y)\right]\tilde{f}(x)$
$\displaystyle=$
$\displaystyle\int\\!\\!dx\,g(x)\tilde{f}(x)-\int\\!\\!dx\\!\\!\int\\!\\!dy\,\tilde{g}(y)K(x,y)\tilde{f}(x)$
$\displaystyle=$
$\displaystyle\int\\!\\!dx\,\tilde{g}(x)\Bigg{[}f(x)-\underbrace{\int\\!\\!dy\,\tilde{f}(y)K(x,y)}_{bf(x)}\Bigg{]}\,,$

where in the second line we used the first relation in Eq. (92), while in the
last line the relabelling $x\leftrightarrow y$ of the integration variables
was carried out, the symmetry of the kernel $K(x,y)$ was exploited, and the
second relation in Eq. (92) was invoked. Thus, one obtains

$\beta=(1-b)\int\\!\\!dx\tilde{f}(x)g(x)\,,$ (95)

where we used that $\tilde{g}(x)f(x)=\tilde{f}(x)g(x)$. Evidently, the
dependence of $\beta$ on the parameter $z$ has been exchanged for the
dependence on the function $g(x)$, which was absent from the original integral
in Eq. (93). In Appendix A, we present a concrete example, where the
equivalence between Eqs. (93) and (95) may be worked out exactly.

Figure 11:  Diagrammatic illustration of the renormalization of $I$. Note that
the diagrams containing the kernel ${\cal K}$ undergo a relabeling of
integration momenta, $k\leftrightarrow\ell$, from the second to the third
line.

At this point it is straightforward to repeat the construction leading to Eq.
(95) for the system of integral equations given by Eqs. (82) and (83), in
conjunction with Eq. (79). To that end, all we need is to establish the
correspondence

$\displaystyle\\{g(x),\,f(x),\,u(x),\,K(x,y)\\}\,$
$\displaystyle\leftrightarrow\,\\{\mathit{L}_{{sg}}(r),\,\mathbb{B}(r),\,k^{2}\Delta^{2}(k),\,\alpha_{s}K(r,k)\\}\,,\,\,$
$\displaystyle\\{\beta,\,z,\,b,\,dx,\,dy\\}\,$
$\displaystyle\leftrightarrow\,\\{2I/3C_{\rm
A},\,Z_{3},\,t,\,d^{4}k,\,d^{4}\ell\\}\,.$ (96)

In particular, with the above identifications and Eq. (71), we find that the
analogue of Eq. (95) is precisely Eq. (90), which is the announced result. The
diagrammatic representation of the steps described in Eq. (94) is shown in
Fig. 11; note that, in doing so, we employ the representation of the BSE given
in Fig. 10, whose main building block was introduced in Fig. 6.

The above construction reveals a striking fact: if the parameter $\omega$ is
set to zero, the cancellation described in Eq. (94) [Fig. 11] is perfect;
thus, even if a non-trivial $\mathbb{B}$ is obtained from Eq. (82), the
renormalized transition amplitude $I$, and therefore, the gluon mass $m$,
vanish. In that sense, the gluon mass emerges due to the mismatch between the
kernels in the equations for $\mathbb{B}$ and $\mathit{L}_{{sg}}(r)$, [see
Eqs. (82) and (83)], produced by the presence of a $t\neq 1$ in the former but
not in the latter. In diagrammatic terms, the non-vanishing of the gluon mass
hinges on the difference between the ${\cal T}$ appearing in Fig. 10 and the
${\cal K}$ in Fig. 8, namely the component ${\cal M}$ [see Fig. 3]. Even
though, on intuitive grounds, the importance of this term seems evident, given
that it carries precisely the information about the emergence of a massless
scalar, the quantitative understanding attained through the above construction
is most valuable.

Quite interestingly, the cancellations described in Eq. (94) are not
coincidental, but are rather driven by an underlying mathematical principle.
In particular, as we discuss in detail in the next section, our SDE-BSE
analysis may be interpreted as a rather subtle application of the so-called
“Fredholm alternatives theorem”.

##  VIII Fredholm alternatives theorem and its evasion

Let us consider the inhomogeneous and homogeneous Fredholm equations of the
second kind, given by [96, 97]

$f_{1}(x)=f_{2}(x)+\lambda\int_{a}^{b}K(x,y)f_{1}(y)dy\,,$ (97)

and

$f_{3}(x)=\lambda\int_{a}^{b}K(x,y)f_{3}(y)dy\,,$ (98)

respectively. Here, $f_{1}$ and $f_{3}$ are the unknown functions, while
$f_{2}$ is a known inhomogeneous term. The kernel, $K(x,y)$, is assumed to be
continuous and square-integrable in the interval $[a,b]$, in which case it is
said to be a Hilbert-Schmidt integral operator.

The Fredholm alternatives theorem imposes restrictions on the existence of
simultaneous solutions of Eqs. (97) and (98). For the purposes of the present
work, it suffices to consider the special case of the theorem when $K(x,y)$ is
real and symmetric in $x\leftrightarrow y$ 555For the generalization of the
theorem for a complex and non-symmetric kernel see, e.g., [96, 97].. Under
these simplifications, the Fredholm alternatives theorem can be stated as
follows [97]:

  1. (a)

If $\lambda$ is not an eigenvalue of $K(x,y)$, i.e., if the homogeneous
equation has only the trivial solution, $f_{3}(x)=0$, then the inhomogeneous
Eq. (97) has a solution for _any_ nonzero $f_{2}(x)$.

  2. (b)

If $\lambda$ is an eigenvalue of $K(x,y)$, such that $f_{3}(x)$ is
nonvanishing, then Eq. (97) has solutions _if and only if_

$\int_{a}^{b}f_{2}(y)f_{3}(y)dy=0\,.$ (99)

In order to connect the Fredholm alternative theorem to the system of
equations satisfied by $\mathit{L}_{{sg}}(r)$ and $\mathbb{B}(r)$, we need to
perform certain transformations that will bring Eqs. (82) and (83) to a form
similar to Eqs. (97) and (98).

We begin by rewriting Eqs. (82) and (83) in spherical coordinates. To this
end, we define the variables

$x:=r^{2}\,,\qquad\qquad y:=k^{2}\,,$ (100)

and the angle $\theta$ between $r$ and $k$, i.e., $r\cdot k=|r||k|c_{\theta}$,
where we use the shorthand notation $c_{\theta}:=\cos\theta$ and
$s_{\theta}:=\sin\theta$. Note that the kernel $K(r,k)$ is a function of $x$,
$y$, and $\theta$, i.e., $K(r,k)\equiv K(x,y,\theta)$.

Moreover, we introduce an ultraviolet momentum cutoff, $\Lambda$, to
regularize potential divergences666For the sake of simplicity, a hard cutoff
has been employed; we have confirmed that exactly the same conclusions are
reached when the calculation is carried out using dimensional regularization.
. In this case, the integral measure takes the form

$\int_{k}:=\frac{1}{(2\pi)^{3}}\int_{0}^{\Lambda^{2}}\\!\\!dy\,y\int_{0}^{\pi}\\!\\!d\theta\,s_{\theta}^{2}\,.$
(101)

Now, the only dependence of the integrands of Eqs. (82) and (83) on the angle
is through $K(x,y,\theta)$. Then, we can define an “angle-integrated kernel”,
${\widehat{K}}(x,y)$, by

${\widehat{K}}(x,y)=\frac{1}{(2\pi)^{3}}\\!\int_{0}^{\pi}\\!\\!d\theta\,s_{\theta}^{2}\,K(x,y,\theta)\,.$
(102)

With the above definitions, Eqs. (82) and (83) are recast as

$\displaystyle\mathbb{B}(x)=$
$\displaystyle\,t^{-1}\alpha_{s}\\!\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\cal
Z}^{2}(y){\widehat{K}}(x,y)\mathbb{B}(y)\,,$
$\displaystyle\mathit{L}_{{sg}}(x)=$
$\displaystyle\,Z_{3}+\alpha_{s}\\!\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\cal
Z}^{2}(y){\widehat{K}}(x,y)\mathit{L}_{{sg}}(y)\,,$ (103)

where we introduced the gluon dressing function, ${\cal Z}(x):=x\Delta(x)$.

At this point, we note that since $K(r,k)$ is symmetric under the exchange of
$r\leftrightarrow k$, then ${\widehat{K}}(x,y)={\widehat{K}}(y,x)$. However,
the complete kernel of Eq. (103), namely ${\cal Z}^{2}(y){\widehat{K}}(x,y)$,
is not symmetric under $x\leftrightarrow y$, due to the factor of ${\cal
Z}^{2}(y)$. Nevertheless, it can be transformed into an equivalent system of
equations with a symmetric kernel.

To this end, we multiply Eq. (103) by ${\cal Z}(x)$, and define

${\widetilde{\mathbb{B}}}(x):={\cal
Z}(x)\mathbb{B}(x)\,,\qquad{\widetilde{L}}_{sg}(x):={\cal
Z}(x)\mathit{L}_{{sg}}(x)\,,\qquad{\widetilde{K}}(x,y):={\cal Z}(x){\cal
Z}(y){\widehat{K}}(x,y)\,.$ (104)

Then, Eq. (103) is equivalent to

$\displaystyle{\widetilde{\mathbb{B}}}(x)=$
$\displaystyle\,t^{-1}\alpha_{s}\\!\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\widetilde{K}}(x,y){\widetilde{\mathbb{B}}}(y)\,,$
$\displaystyle{\widetilde{L}}_{sg}(x)=$ $\displaystyle\,Z_{3}{\cal
Z}(x)+\alpha_{s}\\!\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\widetilde{K}}(x,y){\widetilde{L}}_{sg}(y)\,,$
(105)

whose kernel ${\widetilde{K}}(x,y)$ is symmetric.

Finally, the Eq. (79) for $I$ can be recast as

$I=\frac{3C_{\rm A}Z_{3}}{32\pi^{2}}\\!\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\cal
Z}(y)\,{\widetilde{\mathbb{B}}}(y)\,.$ (106)

The implications of the Fredholm alternatives theorem for Eq. (105) are now
straightforward to establish. Specifically, suppose that $\omega=0$, such that
$t=1$. Then, we have a direct correspondence between Eq. (105) and the Eqs.
(97) and (98) through the identification

$\displaystyle\\{f_{1}(x),\,f_{2}(x),\,f_{3}(x),\,K(x,y)\\}\,$
$\displaystyle\leftrightarrow\,\\{{\widetilde{L}}_{sg}(x),\,Z_{3}{\cal
Z}(x),\,{\widetilde{\mathbb{B}}}(x),\,{\widetilde{K}}(x,y)\\}\,,$
$\displaystyle\\{\lambda,\,a,\,b\\}\,$
$\displaystyle\leftrightarrow\,\\{\alpha_{s},\,0,\,{\Lambda^{2}}\\}\,.$ (107)

Hence, for both $\mathit{L}_{{sg}}(x)$ and $\mathbb{B}(x)$ to be nonzero, the
Fredholm alternatives theorem implies that

$Z_{3}\int_{0}^{\Lambda^{2}}\\!\\!dy\,{\cal
Z}(y)\,{\widetilde{\mathbb{B}}}(y)=0\,.$ (108)

Comparison to Eq. (106) then yields $I=0$, and therefore, due to Eq. (15),
$m^{2}=0$.

We conclude that the generation of a gluon mass through the SM hinges on
$\omega\neq 0$, $t\neq 1$, which leads to the evasion of Fredholm alternatives
theorem, by relaxing the assumption of equal kernels in Eq. (105). Otherwise,
if $\omega=0$, the theorem imposes a vanishing gluon mass, even in the
presence of a nonvanishing $\mathbb{B}$, because, quite remarkably, the
condition of Eq. (99) is precisely the equation for the transition amplitude
$I$.

Let us finally discuss a subtlety related to the applicability of Fredholm’s
theorem in the present context. Specifically, the theorem applies to a system
of _linear_ equations, whereas the kernel $K(r,k)$ depends on
$\mathit{L}_{{sg}}(r)$ [see Fig. 9 and Eq. (115)]. In that sense, even for
$\omega=0$, the equation for $\mathit{L}_{{sg}}(r)$ is nonlinear.
Nevertheless, our main conclusion, namely that for $\omega=0$ the gluon mass
vanishes, persists.

To see this, suppose there exists a solution, $L_{0}(r)$ and
$\mathbb{B}_{0}(r)$, to the full nonlinear system of equations. Then, let
$K_{0}(r,k)$ be the value of $K(r,k)$ obtained by the substitution of
$\mathit{L}_{{sg}}\to L_{0}$ in its expression777As can be seen from Fig. 9,
the $\mathit{L}_{{sg}}$ enters in $K(r,k)$ in such a way that the symmetry
under $r\longleftrightarrow k$ is preserved.. Now, $L_{0}(r)$ and
$\mathbb{B}_{0}(r)$ must also be a solution of

$\displaystyle\mathit{L}_{{sg}}(r)=$
$\displaystyle\,Z_{3}+\alpha_{s}\\!\int_{k}k^{2}\Delta^{2}(k^{2})\,{K}_{0}(r,k)\,\mathit{L}_{{sg}}(k)\,,$
$\displaystyle\mathbb{B}(r)=$
$\displaystyle\,t^{-1}\alpha_{s}\\!\int_{k}k^{2}\Delta^{2}(k^{2})\,{K}_{0}(r,k)\,\mathbb{B}(k)\,,$
(109)

since setting
$\\{\mathit{L}_{{sg}},\,\mathbb{B}\\}\to\\{L_{0},\,\mathbb{B}_{0}\\}$ in the
above equation recovers the original, nonlinear, system. But the Fredholm
alternatives theorem applies to Eq. (109) with $\omega=0$, in which case the
arguments of Sec. VIII lead to $m=0$. Hence, the nonlinear nature of the
equations cannot by itself evade the Fredholm alternatives theorem when
$\omega=0$ and the conclusion of our analysis is unaffected.

##  IX Numerical analysis

In this section we explore the numerical implications of the main results
presented so far, placing particular emphasis on the value of the gluon mass
$m$, and the shape and size of the displacement function ${\mathbb{C}}(r)$. In
doing so, we will use two important results as benchmarks for these
quantities. In particular, we identify as the optimal value for the gluon mass
the inverse of the saturation point of the lattice gluon propagator at the
origin. Evidently, this is a renormalization-point dependent mass scale, and
does not admit a direct physical interpretation888An RGI gluon mass of about
450 MeV may be obtained following the procedure described in [72, 127].; when
the lattice curve has been renormalized such that
$\Delta^{-1}(\mu^{2})=\mu^{2}$ at $\mu=4.3$ GeV, one obtains the value
$m_{{\rm{\scriptscriptstyle lat}}}=354$ MeV [106]. As for the displacement
function, we will compare our results with the curve shown in the left panel
of Fig. 15, which was obtained following the method outlined in Sec. III.4;
for details, see [76].

The main steps of the numerical procedure followed may be summarized as
follows:

(i) Since the four-gluon kernel $K(r,k)$ entering in BSE of Eq. (82) is not
known, we resort to models for it, using its one-gluon exchange approximation,
$K_{\rm oge}(r,k)$, as our point of departure. Specifically, at the beginning
of the numerical procedure, an Ansatz for $K(r,k)$ is constructed, by varying
the parameters of $K_{\rm oge}(r,k)$.

(ii) The $K(r,k)$ from the previous step is fed into the BSE of Eq. (82),
which is solved as an eigenvalue problem, yielding a solution for
$\mathbb{B}(r)$. The corresponding eigenvalue $\alpha_{\star}$ is used in Eq.
(89) to determine the value of $\sigma$, thus fixing the scale of
$\mathbb{B}(r)$. In addition, the value of $\omega$ is obtained directly from
Eq. (87).

(iii) The above solution for $\mathbb{B}(r)$ is substituted into Eq. (90) to
yield the value of $I$. Then, using the two relations in Eq. (84), we get the
value of the gluon mass and the form of the displacement function, and compare
them with the benchmark results.

(vi) The procedure is repeated, using a different model for the kernel
$K(r,k)$.

###  IX.1 One-gluon exchange approximation of the BSE kernel

For practical calculations, the kernel ${\cal
K}_{\rho\sigma\mu\nu}^{mnbc}(-k,k,r,-r)$ of the BSE must be truncated. To this
end, we start with the skeleton expansion of ${\cal
K}_{\rho\sigma\mu\nu}^{mnbc}(-k,k,r,-r)$, shown in Fig. 9, which we write as

${\cal
K}_{\rho\sigma\mu\nu}^{mnbc}(-k,k,r,-r)=\sum_{i=1}^{\infty}(b_{i})_{\rho\sigma\mu\nu}^{mnbc}(-k,k,r,-r)\,.$
(110)

Substituting Eq. (110) into Eq. (81), we can compute the contribution of each
diagram $(b_{i})$ to the Euclidean $K(r,k)$ entering Eq. (82). As it turns
out, the diagram $(b_{2})$, which is simply the tree-level four-gluon vertex,
vanishes when inserted in Eq. (81) [55, 58], i.e.,

$f^{abc}f^{amn}P^{\mu\nu}(r)P^{\rho\sigma}(k)\,(b_{2})^{mnbc}_{\rho\sigma\mu\nu}(-k,k,r,-r)=0\,.$
(111)

Hence, $K(r,k)$ reads

$K(r,k)=K_{\rm oge}(r,k)+K_{\rm ho}(r,k)\,,$ (112)

where $K_{\rm oge}(r,k)$ is the contribution of the one-gluon exchange
diagram, $(b_{1})$, and $K_{\rm ho}(r,k)$ represent the higher-order diagrams
$(b_{i})$ with $i\geq 3$.

Next, we focus on the one-gluon exchange contribution, $K_{\rm oge}(r,k)$.
From diagram $(b_{1})$ of Fig. 9 and Eq. (81), one obtains in the Landau gauge

$K_{\rm oge}(r,k)=\frac{2\pi C_{\rm
A}}{3}\Delta(k-r)\left[\overline{\Gamma}_{\mu\rho\sigma}(r,-k,k-r)\overline{\Gamma}^{\mu\rho\sigma}(-r,k,r-k)\right]_{{\rm{\scriptscriptstyle
E}}}\,,$ (113)

where

$\overline{\Gamma}_{\mu\rho\sigma}(r,-k,k-r):=P^{\mu^{\prime}}_{\mu}(r)P^{\rho^{\prime}}_{\rho}(k)P^{\sigma^{\prime}}_{\sigma}(k-r){\rm{I}}\\!\Gamma_{\mu^{\prime}\rho^{\prime}\sigma^{\prime}}(r,-k,k-r)\,,$
(114)

is the transversely projected three-gluon vertex [128, 129, 130, 73, 74, 105,
110, 116, 53].

In principle, the $\overline{\Gamma}_{\mu\rho\sigma}$ appearing in Eq. (113)
consists of four independent tensor structures, each accompanied by a form
factor that depends on three kinematic variables. However, as has been shown
in numerous recent works [128, 129, 130, 73, 74, 105, 131, 110, 116, 53], the
classical tensor structure of $\overline{\Gamma}_{\mu\rho\sigma}$ is dominant,
and the associated form factor $\mathit{L}_{{sg}}$ may be accurately described
as a function of a single special kinematic variable, denoted by $s^{2}$,
namely

$\overline{\Gamma}^{\mu\rho\sigma}(r,-k,k-r)=\overline{\Gamma}_{\\!0}^{\mu\rho\sigma}(r,-k,k-r)\mathit{L}_{{sg}}(s)\,,\qquad
s^{2}:=\frac{1}{2}[r^{2}+k^{2}+(k-r)^{2}]\,,$ (115)

where $\overline{\Gamma}_{\\!0}^{\mu\rho\sigma}$ is the tree-level value of
$\overline{\Gamma}^{\mu\rho\sigma}$, obtained by substituting the
${\rm{I}}\\!\Gamma^{\mu^{\prime}\rho^{\prime}\sigma^{\prime}}$ in Eq. (114) by
the ${\rm{I}}\\!\Gamma_{\\!0}^{\mu^{\prime}\rho^{\prime}\sigma^{\prime}}$ of
Eq. (5).

Then, using Eq. (115), employing spherical coordinates, and introducing the
variables defined in Eq. (100), one obtains

$K_{\rm oge}(x,y,\theta)=K_{\rm{kin}}(x,y,\theta){\cal R}_{\rm oge}(u,s)\,,$
(116)

where $u^{2}:=(k-r)^{2}=x+y-2c_{\theta}\sqrt{xy}$,
$s^{2}=x+y-c_{\theta}\sqrt{xy}$,

$K_{\rm{kin}}(x,y,\theta):=\frac{8\pi C_{\rm
A}}{3u^{2}\sqrt{xy}}\left\\{c_{\theta}s_{\theta}^{2}\left[(c_{\theta}^{2}+8)xy-6c_{\theta}\sqrt{xy}(x+y)+3(x^{2}+y^{2})\right]\right\\}\,,$
(117)

and

${\cal R}_{\rm oge}(u,s):={\Delta(u)\mathit{L}_{{sg}}^{2}(s)}\,.$ (118)

The ingredients entering in the $K_{\rm oge}(x,y,\theta)$ of Eq. (116) are all
accurately known from lattice simulations. In particular, we use for
$\Delta(u)$ and $\mathit{L}_{{sg}}(s)$ the fits to the lattice results of
[106], given by Eqs. (C.11) and (C.12) of [59], respectively; the
corresponding curves are shown in Fig. 12. Note that these ingredients are
renormalized in the so-called “asymmetric MOM scheme” [102, 103, 105, 106,
116], defined by the renormalization condition

$\Delta^{-1}(\mu)=\mu^{2}\,,\qquad\mathit{L}_{{sg}}(\mu)=1\,,$ (119)

with $\mu=4.3$ GeV denoting the renormalization point. For this
renormalization scheme and point, the corresponding strong charge takes the
value $\alpha_{s}=0.27$ [103], which we adopt from now on.

Figure 12: Lattice data (points) for the gluon propagator from [106] (left)
and the three-gluon vertex form factor, $\mathit{L}_{{sg}}(s)$, from [105]
(right). The blue continuous curves represent the corresponding fits given by
Eqs. (C11) and (C12) of [59].

###  IX.2 Asymptotic behavior of $\mathbb{B}(r)$

For a given kernel, $K(r,k)$, the eigenvalue problem of Eq. (86) can be solved
with standard procedures, such as _Nyström_ ’s method [132]. This determines
the eigenvalue, $\alpha_{\star}$, and an eigenfunction, $\mathbb{B}_{0}(r)$,
arbitrarily normalized such that its global maximum is $1$. In order to fix
the physical scale of $\mathbb{B}(r)$ and compute the gluon mass, we then need
to compute $\omega$ and $I$ through Eqs. (78) and (90). This task requires an
analysis of the ultraviolet behavior of the integrals in the latter equations,
which we describe in detail below. For simplicity, we will illustrate the
procedure with the case of the one-gluon exchange kernel, $K_{\rm oge}$, but
it applies equally to the kernels constructed through Eq. (126).

Let us begin by writing Eqs. (78) and (90) in spherical coordinates,

$\displaystyle\omega=$
$\displaystyle\,\lambda\int_{0}^{\infty}\\!\\!\\!dy\,y^{2}{\Delta}^{2}(y)\,\mathbb{B}^{2}(y)\,,$
(120) $\displaystyle I=$
$\displaystyle\,\lambda\omega\int_{0}^{\infty}\\!\\!\\!dy\,y^{2}{\Delta}^{2}(y)\,\mathit{L}_{{sg}}(y)\,\mathbb{B}(y)\,,$
(121)

where $\lambda:=3C_{\rm A}/(32\pi^{2})$.

In addition, let us recall that at large momenta $\Delta(y)$ and
$\mathit{L}_{{sg}}(y)$ approach their one-loop resumed perturbative behaviors
[133, 60, 134, 135, 136, 73]

$\Delta(y)\sim y^{-1}L_{{\rm{\scriptscriptstyle
UV}}}^{-13/22}(y)\,,\qquad\mathit{L}_{{sg}}(y)\sim L_{{\rm{\scriptscriptstyle
UV}}}^{17/44}(y)\,,\qquad L_{{\rm{\scriptscriptstyle
UV}}}(y):=c\ln(y/\Lambda^{2}_{{\rm{\scriptscriptstyle MOM}}})\,,$ (122)

where $c=1/\ln(\mu^{2}/\Lambda^{2}_{{\rm{\scriptscriptstyle MOM}}})$, and
$\Lambda_{{\rm{\scriptscriptstyle MOM}}}$ is the (quenched) QCD mass-scale in
the MOM scheme [137, 138, 139, 140]; we employ the one-loop result
$\Lambda^{2}_{{\rm{\scriptscriptstyle MOM}}}=\mu^{2}\exp[-12\pi/(11C_{\rm
A}\alpha_{s})]$, for which $\Lambda=520$ MeV and $c=0.236$. Note that the fits
of Eqs. (C.11) and (C.12) of [59] for $\Delta(r)$ and $\mathit{L}_{{sg}}(r)$
reproduce Eq. (122) by construction.

Figure 13: Absolute value of the eigenfunction $\mathbb{B}_{0}(r)$ obtained
with the one-gluon exchange kernel, $K_{\rm oge}(r,k)$, normalized such that
it has $1$ as global maximum (blue continuous). The orange dot-dashed line
corresponds to the asymptotic form of Eq. (123), and accurately approximates
$\mathbb{B}_{0}(r)$ for $r^{2}>25$ GeV2.

Now, the asymptotic behavior of $\mathbb{B}_{0}(r)$ is not known analytically.
Nevertheless, the numerical solutions for $\mathbb{B}_{0}(r)$ are found to be
accurately approximated by the asymptotic form

$\lim_{r\to\infty}\mathbb{B}_{0}(r)=\frac{a}{r^{2}L_{{\rm{\scriptscriptstyle
UV}}}^{\kappa}(r^{2})}\,,$ (123)

for $r>5$ GeV, where the constants $a$ and $\kappa$ can be determined by
fitting the $\mathbb{B}_{0}(r)$ data in this range. In Fig. 13 this is shown
explicitly for the case when $K(r,k)=K_{\rm oge}(r,k)$, with the blue
continuous curve showing the full solution $\mathbb{B}_{0}(r)$, whereas the
orange dot-dashed represents Eq. (123) with $a=3.17$ and $\kappa=0.958$. The
corresponding eigenvalue is $\alpha_{\star}=0.685$.

Next, to get $\omega_{0}$, we set $\mathbb{B}\to\mathbb{B}_{0}$ into Eq.
(120). With the aid of Eqs. (122) and (123), it is easy to establish that the
integral converges rapidly, and its numerical evaluation is rather
straightforward, yielding $\omega_{0}=0.389$.

In contrast, the expression for $I$ given in Eq. (121) is found to converge
slowly, requiring a more careful treatment.

To analyze the ultraviolet convergence of Eq. (121), we first separate the
integration interval in two parts, by means of a large momentum scale,
$M^{2}$, i.e.,

$I=\underbrace{\lambda\omega\int_{0}^{M^{2}}\\!\\!\\!dy\,y^{2}{\Delta}^{2}(y)\,\mathit{L}_{{sg}}(y)\,\mathbb{B}(y)}_{I_{{\rm{\scriptscriptstyle
IR}}}}+\underbrace{\lambda\omega\int_{M^{2}}^{\infty}\\!\\!\\!dy\,y^{2}{\Delta}^{2}(y)\,\mathit{L}_{{sg}}(y)\,\mathbb{B}(y)}_{I_{{\rm{\scriptscriptstyle
UV}}}}\,.$ (124)

For $M$ sufficiently large, all of the inputs in $I_{{\rm{\scriptscriptstyle
UV}}}$ can be substituted by their asymptotic forms, namely Eqs. (122) and
(123). The resulting integral can then be done analytically,

$\displaystyle I_{{\rm{\scriptscriptstyle UV}}}\sim$ $\displaystyle\,\sigma
a\lambda\omega\int_{M^{2}}^{\infty}\\!\\!\\!dy\,\frac{1}{yL^{\kappa+35/44}_{{\rm{\scriptscriptstyle
UV}}}(y)}=\frac{\sigma
a\lambda\omega}{c^{\kappa+35/44}(\kappa-9/44)\ln^{\kappa-9/44}(M^{2}/\Lambda^{2}_{{\rm{\scriptscriptstyle
MOM}}})}\,,$ (125)

provided that $\kappa>9/44=0.204$, which we find to be satisfied by our
numerical solutions for $\mathbb{B}(r)$.

However, the convergence of Eq. (125) is slow, and direct numerical
integration of Eq. (121) becomes unreliable. Instead, in what follows we use
Eq. (124) with $M=5$ GeV, and compute $I$ in two steps: first,
$I_{{\rm{\scriptscriptstyle IR}}}$ is calculated by direct numerical
integration of the interpolated data for $\mathbb{B}(r)$; then,
$I_{{\rm{\scriptscriptstyle UV}}}$ is approximated by Eq. (125), and evaluated
with the previously determined values of $a$, $\kappa$ and $\sigma$. Finally,
this procedure yields $I=\pm 0.689$ GeV.

Then, using Eq. (84), we obtain $m=1.27$ GeV, and the resulting displacement
function ${\mathbb{C}}(r)$ is shown as a yellow dashed curve in Fig. 15.
Evidently, both findings are considerably larger than the results obtained
from the lattice analysis of [106, 76] for $\\{m_{{\rm{\scriptscriptstyle
lat}}},{\mathbb{C}}_{{\rm{\scriptscriptstyle lat}}}(r)\\}$. This, in turn,
motivates the study of the potential impact that the form of the kernel may
have on these quantities.

###  IX.3 Varying the BSE kernel

In this subsection, we modify the form of the BSE kernel for the purpose of
reaching a better agreement with the lattice results.

It is clear that, in order to determine the effect of the higher-order
diagrams comprising $K_{\rm ho}(r,k)$, one has to carry out extensive
calculations. Instead, in the present analysis we simply treat the full
$K(r,k)$ as a deformation of the one-gluon exchange $K_{\rm oge}(r,k)$,
obtained by varying the function $\Delta(u)$ entering in Eq. (116). We
emphasize that this is the only place where this function is varied;
everywhere else, $\Delta$ retains its standard form.

Specifically, in the kernel of the BSE for $\mathbb{B}(r)$, see Eq. (101), we
set

$K(x,y,\theta)\to K_{\rm eff}(x,y,\theta)=K_{\rm{kin}}(x,y,\theta){\cal
R}_{\rm eff}(u,s)\,,$ (126)

where the function ${\cal R}_{\rm eff}(u,s)$ is given by

${\cal R}_{\rm eff}(u,s)=\Delta_{\rm eff}(u)\mathit{L}_{{sg}}(s)\,,$ (127)

with $\Delta_{\rm eff}(u)$ parametrized as

$\Delta_{\rm
eff}(u)=\Delta(u)\times\left[1+\frac{c_{0}u^{2}}{1+c_{1}u^{2}+c_{2}u^{4}}\right]\,.$
(128)

Note that $K_{\rm oge}$ is recovered by setting $c_{0}=0$ in Eq. (128).
Moreover, at large momenta, $\Delta_{\rm eff}(u)$ reduces to $\Delta(u)$, such
that the $K_{\rm eff}$ reduce to $K_{\rm oge}$.

We next vary the $c_{i}$ in Eq. (128) within certain intervals, and consider
the resulting behaviour of $\\{m,\,{\mathbb{C}}\\}$. Upon systematic
inspection of the parameter space, we conclude that the best values for
$c_{1}$ and $c_{2}$ are $c_{1}=0.00667$ GeV-2 and $c_{2}=0.0486$ GeV-4; we
therefore keep $c_{1}$ and $c_{2}$ fixed at these values, and only vary
$c_{0}$.

Figure 14:  Top left: A family of variations of $\Delta_{\rm eff}(u)$,
obtained by varying the parameter $c_{0}$ in Eq. (128). The inset shows the
corresponding dressing functions, ${\cal Z}_{\rm{eff}}(q):=q^{2}\Delta_{\rm
eff}(q)$. Top right: Diagonal ${\widehat{K}}_{\rm{eff}}(k^{2},k^{2})$ of the
angular-integrated kernel ${\widehat{K}}_{\rm{eff}}(r^{2},k^{2})$, defined by
Eq. (102), corresponding to the family of $\Delta_{\rm eff}(u)$ on the left
panel. Bottom: General kinematics ${\widehat{K}}(r^{2},k^{2})$ for $K\to
K_{\rm oge}$ (left), $K\to K_{\rm{eff}}$ (right). Note that the maximum of
$K_{\rm{eff}}$ is enhanced by about $30\%$ with respect to that of
${\widehat{K}}_{\rm{oge}}$.

On the top left panel of Fig. 14 we show the family of $\Delta_{\rm eff}(r)$
obtained by varying $c_{0}$ within the interval $[0,0.503]$ GeV-2, while on
the top right panel the diagonals ($r^{2}=k^{2}$) of the resulting
${\widehat{K}}_{\rm{eff}}(r^{2},k^{2})$ are displayed. As we see, enhancing
$\Delta_{\rm eff}(u)$ leads to an enhancement of
${\widehat{K}}_{\rm{eff}}(r^{2},k^{2})$. In the bottom panel one may see how
the entire shape of the angular-integrated kernel gets modified when $c_{0}$
varies from $c_{0}=0$, corresponding to the
${\widehat{K}}_{\rm{oge}}(r^{2},k^{2})$, to the case $c_{0}=0.503$ GeV-2.
Note, in particular, that the maximum of
${\widehat{K}}_{\rm{eff}}(r^{2},k^{2})$ is enhanced by about $30\%$ with
respect to that of ${\widehat{K}}_{\rm{oge}}$.

Next, we compute the values of $m$ and ${\mathbb{C}}(r)$ corresponding to the
family of $\Delta_{\rm eff}(u)$ of Fig. 14. The resulting displacement
amplitudes are shown on the left panel of Fig. 15, with the corresponding
values of $m$ given in the legend.

Note that, for $c_{0}=0.212$ GeV-2 (red continuous curve on the panel of Fig.
15), the resulting ${\mathbb{C}}(r)$ is statistically consistent with
${\mathbb{C}}_{{\rm{\scriptscriptstyle lat}}}(r)$; however, for this value of
$c_{0}$, the mass turns out to be $m=687$ MeV, nearly twice the value of
$m_{{\rm{\scriptscriptstyle lat}}}$. Variation of the remaining parameters of
Eq. (128) show the same pattern: improving the agreement between
${\mathbb{C}}(r)$ and ${\mathbb{C}}_{{\rm{\scriptscriptstyle lat}}}(r)$ causes
$m$ to exceed $m_{{\rm{\scriptscriptstyle lat}}}$.

Figure 15:  Left: Displacement amplitude, ${\mathbb{C}}(r)$, for different
values of the parameter $c_{0}$ (see legends of Fig. 14) of Eq. (128). The
lattice-derived results of [76] are shown as points for comparison. Note that
the red dashed line lies entirely within the errors of
${\mathbb{C}}_{{\rm{\scriptscriptstyle lat}}}(r)$. Right: Gluon propagator,
$\Delta(q)$, for different masses, assuming the functional form of Eq. (129),
compared to the lattice data of [106].

Lastly, we illustrate the effect of the masses obtained in the Fig. 15 on the
gluon propagator. Evidently, determining $\Delta(q)$ self-consistently
requires solving the coupled system of equations comprised by the gluon
propagator and three-gluon vertex SDEs, which is beyond the scope of the
present study. Instead, we assume that $\Delta(q)$ is modified only near the
origin by the value of $m$, recovering the lattice results for large $q$. A
functional form that implements this assumption is

$\Delta(q)=\Delta_{{\rm{\scriptscriptstyle
fit}}}(q)-[\Delta_{{\rm{\scriptscriptstyle fit}}}(0)-1/m^{2}]\exp(-q/\nu)\,,$
(129)

where $\nu=0.05$ GeV, and $\Delta_{{\rm{\scriptscriptstyle fit}}}(q)$ is the
fit to the lattice data of [106] given by Eq. (C.11) of [59]. The resulting
curves are shown on the right panel of Fig. 15.

##  X Discussion and conclusions

In the present work we have addressed important open issues related to the
implementation of the SM in a non-Abelian context. In particular, the question
of how the scale of the gluon mass may be fixed has been solved, at least in
principle, by introducing a nonlinear component into the BSE that controls the
bound state amplitude, $\mathbb{B}(r)$, of the massless pole formation. In
addition, we have put forth a novel approach for dealing with the problem of
the multiplicative renormalizability of the gluon mass equation, thus reaching
a manifestly finite answer for the gluon mass. The success of this
renormalization procedure hinges on the implementation of a key cancellation,
which proceeds through the judicious combination of the SDE for the three-
gluon vertex and the BSE for $\mathbb{B}(r)$. It turns out that the
mathematical origin of this cancellation may be traced back to Fredholm’s
alternatives theorem, operating at the level of the SDE-BSE system that we
consider. The nonlinear character of the BSE is crucial for evading this
theorem, which would otherwise enforce the strict vanishing of the gluon mass.

For the purpose of simplifying this rather technical subject, in our
presentation we have considered only the minimum of diagrams required for
exposing the crucial cancellations [see Sec. VII] and their connection to
Fredholm’s theorem. It is, however, of the utmost importance to emphasize
that, as we have explicitly confirmed, the omitted contributions undergo
themselves completely analogous cancellations, which proceed in exactly the
same way, and for precisely the same reason. Actually, the construction
related to the renormalization may be repeated, involving the ghost-gluon and
four-gluon vertices, their SDE-BSE systems, and the corresponding vertex
renormalization constants associated with them. In fact, the only
contributions that survive, and may induce modifications to the results
presented here originate from the non-linear components of the additional BSEs
that control the formation of poles in the ghost-gluon and four-gluon
vertices. The detailed elaboration of these issues will be presented in a
forthcoming communication.

The numerical treatment presented in Sec. IX focused on two principal
quantities, namely the value of the gluon mass, $m$, and the shape and size of
the displacement function, ${\mathbb{C}}(r)$. The procedure followed consists
in producing a sequence of models for $K(r,k)$, in an attempt to optimize the
resulting set $\\{m,\,{\mathbb{C}}(r)\\}$ with respect to the lattice-based
benchmarks for these quantities. The analysis reveals a preference towards
gluon masses that are larger than what the saturation point of the lattice
gluon propagator would indicate. In particular, the red continuous curve in
Fig. 15, which shows the best coincidence with the ${\mathbb{C}}(r)$ of [76],
is obtained for $m=687$ MeV, which is to be contrasted to the
$m_{{\rm{\scriptscriptstyle lat}}}=354$ MeV [106]. As we discuss below, there
are at least three reasons that could account for this relative tension
between $m$ and ${\mathbb{C}}(r)$.

First, for our numerical analysis we have essentially resorted to models for
the kernel, using the one-gluon exchange version as our reference. Evidently,
the next important task in this context would be to actually compute the
modifications induced to $K(r,k)$ by the inclusions of the one-loop dressed
diagrams ($b_{3}$)-($b_{6}$) appearing in the skeleton expansion of Fig. 9.
The result of such a computation is likely to modify the shape of
$\mathbb{B}(r)$, in ways not captured by the modeling described by Eqs. (126)
and (128), thus improving the matching between $m$ and ${\mathbb{C}}(r)$.
Calculations in this direction are already underway; we hope to report on the
results of this study in the near future.

Second, the indirect determination of ${\mathbb{C}}(r)$ presented in [76]
involves a particular partial derivative of the ghost-gluon kernel,
$H_{\nu\mu}(p,q,r)$, appearing in the STI of Eq. (39). This derivative has
never been simulated on the lattice; in the analysis of [76] it was estimated
with the aid of the “one-loop dressed” SDE that controls the evolution of
$H_{\nu\mu}(p,q,r)$. Preliminary analysis indicates that a relatively small
change in the functional form of this derivative would suffice for bringing
$m$ and ${\mathbb{C}}(r)$ to a closer agreement. It is therefore important to
examine the modifications induced to this quantity when further contributions
to the SDE are included; in particular, such a study requires the inclusion
into the aforementioned SDE of the ghost-ghost-gluon-gluon vertex, first
explored in [141].

Third, as mentioned above, the inclusion of the remaining fundamental vertices
of the theory, namely the ghost-gluon and four-gluon vertices, is likely to
modify the expression for $I$ given in Eq. (90). In particular, the non-
linearities in the BSEs that control the formation of the composite scalar
through the fusion of a ghost-antighost pair or three gluons will provide
additional finite contributions to the $I$, thus affecting the value of the
gluon mass. The inclusion of the ghost-gluon vertex and the associated pole
dynamics is well within our grasp, and various of its main aspects have
already been addressed in the literature [75, 58, 59, 49]. On the other hand,
the pole content and structure of the four-gluon vertex are largely
unexplored, and quantitative information is rather difficult to obtain [142,
143, 144, 73, 145, 146]; nonetheless, the understanding acquired from the
present analysis is expected to facilitate future attempts in this direction.

##  XI Acknowledgments

The work of M. N. F. is supported by the National Natural Science Foundation
of China (grant no. 12135007). The research of J. P. is supported by the
Spanish MICINN grant PID2020-113334GB-I00, the Generalitat Valenciana grant
CIPROM/2022/66, and in part by the EMMI visiting grant of the ExtreMe Matter
Institute EMMI at the GSI Helmholtzzentrum für Schwerionenforschung GmbH,
Darmstadt, Germany. J. P. thanks J. M. Pawlowski for his kind hospitality at
the Institute for Theoretical Physics, University of Heidelberg, where part of
this work was completed, and for numerous stimulating discussions.

##  Appendix A A toy model

In this Appendix we present a concrete example of the mathematical
construction described in Sec.VII, leading from Eq. (92) to Eq. (95).

In order to make the analogy with the equations for $\mathit{L}_{{sg}}(r)$,
$\mathbb{B}(r)$, and $I$ as realistic as possible, we use an integration
interval from $0$ to $\Lambda^{2}$, with $\Lambda$ an ultraviolet cutoff,
which will eventually be taken to infinity.

For simplicity, we consider the case of a separable kernel, $K(x,y)=A(x)B(y)$,
for which Eq. (92) can be solved formally. Then, the symmetry of the kernel
under the exchange of $x$ and $y$ implies $B(y)=A(y)$. Moreover, we set
$u(x)=1$. In this case, Eq. (92) simplifies to

$g(x)=z+A(x)\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,g(y)A(y)\,,\qquad
f(x)=b^{-1}A(x)\\!\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,f(y)A(y)\,.$ (130)

To fix the ideas, let us first ignore convergence issues of the integrals in
Eq. (130) as $\Lambda\to\infty$. Then, the formal solution of Eq. (130) for
$g(x)$ and $f(x)$ is read off immediately,

$g(x)=z+cA(x)\,,\qquad\qquad f(x)=dA(x)\,,$ (131)

where $c$ and $d$ are so-far unknown constants, defined as

$c=\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,g(y)A(y)\,,\qquad\qquad
d=b^{-1}\\!\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,f(y)A(y)\,.$ (132)

Next, the constant $c$ can be determined by substituting Eq. (131) into the
first of Eq. (132) to eliminate $g(y)$. This yields a linear equation for $c$,
whose solution reads

$c=z\left[1-\int_{0}^{\Lambda^{2}}\\!\\!\\!dw\,A^{2}(w)\right]^{-1}\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,A(y)\,.$
(133)

On the other hand, substituting Eq. (131) into the second of Eq. (132) leads
to

$d=db^{-1}\\!\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,A^{2}(y)\,,$ (134)

which leaves $d$ undetermined. That is to be expected, since the equation for
$f(x)$ is homogeneous; its solutions can only be determined up to a
multiplicative constant. Instead, Eq. (134) determines $b$, which plays the
role of eigenvalue. Specifically,

$b=\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,A^{2}(y)\,.$ (135)

Now we consider the two forms of $\beta$, given by Eqs. (93) and (95). Using
Eq. (131), the Eqs. (93) and (95) reduce to

$\beta=z\,d\\!\int_{0}^{\Lambda^{2}}\\!\\!dx\,A(x)\,,\qquad\qquad\beta=d\,(1-b)\\!\int_{0}^{\Lambda^{2}}\\!\\!dx\,A(x)[z+cA(x)]\,.$
(136)

Then, using Eqs. (133) and (135) for $c$ and $b$, one finds by simple algebra
that the two expressions yield the same result.

To make this example more concrete, let $A(x)=e^{-x}$. With this kernel, all
of the above integrals converge in the limit $\Lambda\to\infty$, which we can
take directly. Then, we can set $z=1$, such that Eq. (92) reduces to

$g(x)=1+\int_{0}^{\infty}\\!\\!\\!dy\,g(y)e^{-x-y}\,,\qquad\qquad
f(x)=b^{-1}\\!\int_{0}^{\infty}\\!\\!\\!dy\,f(y)e^{-x-y}\,.$ (137)

The solution is found from Eqs. (131), (133), and (135), which yield $b=1/2$,
$c=2$, and

$g(x)=1+2e^{-x}\,,\qquad\qquad f(x)=de^{-x}\,,$ (138)

with $d$ an undetermined constant. These results are readily verified to
satisfy Eq. (137). Finally, both forms of Eq. (136) yield $\beta=d$.

Now, the actual equations for $\mathit{L}_{{sg}}(r)$ and $\mathbb{B}(r)$
contain ultraviolet divergences that are cured by renormalization. To emulate
this situation within our toy model, let us choose

$A(x)=\frac{a}{x+h^{2}}\,,$ (139)

with $a$ and $h$ constant. With this kernel, no infrared divergences arise,
but the treatment of the $\Lambda\to\infty$ limit requires additional care.

The homogeneous equation for $f(x)$ does not present a problem; its solution
is still given by Eq. (138), with the ultraviolet finite eigenvalue

$b=\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,A^{2}(y)=\frac{a^{2}\Lambda^{2}}{h^{2}(\Lambda^{2}+h^{2})}\,,\qquad\qquad\lim_{\Lambda\to\infty}b=\frac{a^{2}}{h^{2}}\,.$
(140)

However, if $z$ were cutoff-independent, say $z=1$, the inhomogeneous equation
for $g(x)$ in Eq. (130) would lead to a divergent result. Specifically, from
Eq. (133),

$c=\frac{ah^{2}(\Lambda^{2}+h^{2})}{(h^{2}-a^{2})\Lambda^{2}+h^{4}}\ln\left(\frac{\Lambda^{2}+h^{2}}{h^{2}}\right)\,,$
(141)

such that

$\lim_{\Lambda\to\infty}c=\frac{ah^{2}}{h^{2}-a^{2}}\ln\left(\frac{\Lambda^{2}}{h^{2}}\right)\,,$
(142)

i.e., $c$, and therefore $g(x)$, diverge logarithmically, in close analogy to
the unrenormalized $\mathit{L}_{{sg}}(x)$.

Let us then choose $z$ in Eq. (130) to be cutoff-dependent, in such a way that
it cancels the divergence of the integral in the equation for $g(x)$. To this
end, we impose a “renormalization condition”, $g(\mu^{2})=1$, at some
renormalization point, $\mu^{2}$. This condition determines $z$ as

$z=1-A(\mu^{2})\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,g(y)A(y)\,.$ (143)

Then, the first of Eq. (92) becomes

$g(x)=1+\left[A(x)-A(\mu^{2})\right]\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,g(y)A(y)\,,$
(144)

while the equation for $f(x)$ is unchanged.

The solution of the system is reached through the same steps used to obtain
Eqs. (131) and (132). Specifically, we find

$g(x)=1+c_{{\rm{\scriptscriptstyle R}}}[A(x)-A(\mu^{2})]\,,\qquad\qquad
f(x)=d_{{\rm{\scriptscriptstyle R}}}A(x)\,,$ (145)

with

$\displaystyle c_{{\rm{\scriptscriptstyle R}}}=$
$\displaystyle\,\left\\{1-\int_{0}^{\Lambda^{2}}\\!\\!\\!dw\,A(w)\left[A(w)-A(\mu^{2})\right]\right\\}^{-1}\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,A(y)$
$\displaystyle=$
$\displaystyle\,\frac{ah^{2}\left(\Lambda^{2}+h^{2}\right)\left(\mu^{2}+h^{2}\right)\ln\left(1+\Lambda^{2}/h^{2}\right)}{\left(\mu^{2}+h^{2}\right)\left[\Lambda^{2}(h^{2}-a^{2})+h^{4}\right]+a^{2}h^{2}\left(\Lambda^{2}+h^{2}\right)\ln\left(1+\Lambda^{2}/h^{2}\right)}\,,$
(146)

$d_{{\rm{\scriptscriptstyle R}}}$ arbitrary, and the eigenvalue $b$ given by
Eq. (140).

Now, the constant $c_{{\rm{\scriptscriptstyle R}}}$, and hence $g(x)$, is
finite. In fact, for $\Lambda\to\infty$, $c_{{\rm{\scriptscriptstyle R}}}$
takes the simple form

$\lim_{\Lambda\to\infty}c_{{\rm{\scriptscriptstyle
R}}}=\frac{\mu^{2}+h^{2}}{a}\,.$ (147)

Next, we consider the two expressions for $\beta$ in Eq. (136). To calculate
the first form, we first need to know the explicit value of $z$; this is
achieved by substituting the solution given in Eqs. (145) and (146) for $g(x)$
into Eq. (143), which furnishes

$z=\frac{\left(\mu^{2}+h^{2}\right)\left[\Lambda^{2}(h^{2}-a^{2})+h^{4}\right]}{\left(\mu^{2}+h^{2}\right)\left[\Lambda^{2}(h^{2}-a^{2})+h^{4}\right]+a^{2}h^{2}\left(\Lambda^{2}+h^{2}\right)\ln\left(1+\Lambda^{2}/h^{2}\right)}\,.$
(148)

Then, with a little algebra, both expressions in Eq. (136) can be shown to
yield the same result, namely

$\beta=\frac{ad\left(\mu^{2}+h^{2}\right)\left[\Lambda^{2}(h^{2}-a^{2})+h^{4}\right]\ln\left(1+\Lambda^{2}/h^{2}\right)}{\left(\mu^{2}+h^{2}\right)\left[\Lambda^{2}(h^{2}-a^{2})+h^{4}\right]+a^{2}h^{2}\left(\Lambda^{2}+h^{2}\right)\ln\left(1+\Lambda^{2}/h^{2}\right)}\,,$
(149)

which for $\Lambda\to\infty$ has the finite value

$\lim_{\Lambda\to\infty}\beta=\frac{d(\mu^{2}+h^{2})(h^{2}-a^{2})}{ah^{2}}\,.$
(150)

We conclude this exercise by considering the scale-fixing of the homogeneous
equation; this can be achieved by imposing an additional relation between $b$
and $f$.

To retain a close analogy with the situation encountered in the equations for
$\mathbb{B}(r)$, we recall that, in our toy model, $b$ plays the role of the
$t$ appearing in Eq. (76). Then, by analogy to Eq. (71), we impose that $b$
must satisfy

$b=1-\int_{0}^{\Lambda^{2}}\\!\\!\\!dy\,f^{2}(y)\,,$ (151)

with the integral on the r.h.s. emulating the $\omega$ of Eq. (71). Note that
since Eq. (151) is quadratic in $f(y)$, Eqs. (139) and (145) yield an
ultraviolet finite $b$, namely

$b=1-\frac{d_{{\rm{\scriptscriptstyle
R}}}^{2}a^{2}\Lambda^{2}}{h^{2}(\Lambda^{2}+h^{2})}\,,\qquad\qquad\lim_{\Lambda\to\infty}b=1-\frac{d_{{\rm{\scriptscriptstyle
R}}}^{2}a^{2}}{h^{2}}\,.$ (152)

Now, Eqs. (140) and (152) must yield the same value of $b$. Hence, equating
them fixes the scale of $d_{\rm R}$ to

$d_{{\rm{\scriptscriptstyle
R}}}=\pm\sqrt{\frac{h^{2}(\Lambda^{2}+h^{2})}{a^{2}\Lambda^{2}}-1}\,,\qquad\qquad\lim_{\Lambda\to\infty}d_{{\rm{\scriptscriptstyle
R}}}=\pm\sqrt{\frac{h^{2}}{a^{2}}-1}\,.$ (153)

Since Eq. (151) is quadratic in $f(y)$, the sign of
$d_{{\rm{\scriptscriptstyle R}}}$ is left undetermined, similarly to the sign
of $\sigma$ in Eq. (89).

## References

  * Yang and Mills [1954] C.-N. Yang and R. L. Mills, Phys. Rev. 96, 191 (1954).
  * Gross and Wilczek [1973] D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343 (1973).
  * Politzer [1973] H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973).
  * Marciano and Pagels [1978] W. J. Marciano and H. Pagels, Phys. Rept. 36, 137 (1978).
  * Faddeev and Popov [1967] L. D. Faddeev and V. N. Popov, Phys. Lett. B 25, 29 (1967).
  * Becchi  _et al._ [1975] C. Becchi, A. Rouet, and R. Stora, Commun. Math. Phys. 42, 127 (1975).
  * Tyutin [1975] I. V. Tyutin, LEBEDEV-75-39 (1975).
  * Becchi  _et al._ [1976] C. Becchi, A. Rouet, and R. Stora, Annals Phys. 98, 287 (1976).
  * ’t Hooft and Veltman [1972] G. ’t Hooft and M. J. G. Veltman, Nucl. Phys. B 44, 189 (1972).
  * Bollini and Giambiagi [1972] C. G. Bollini and J. J. Giambiagi, Nuovo Cim. B 12, 20 (1972).
  * Cornwall [1982] J. M. Cornwall, Phys. Rev. D 26, 1453 (1982).
  * Schwinger [1962a] J. S. Schwinger, Phys. Rev. 125, 397 (1962a).
  * Schwinger [1962b] J. S. Schwinger, Phys. Rev. 128, 2425 (1962b).
  * Eichten and Feinberg [1974] E. Eichten and F. Feinberg, Phys. Rev. D 10, 3254 (1974).
  * Smit [1974] J. Smit, Phys. Rev. D 10, 2473 (1974).
  * Parisi and Petronzio [1980] G. Parisi and R. Petronzio, Phys. Lett. B 94, 51 (1980).
  * Bernard [1982] C. W. Bernard, Phys. Lett. B 108, 431 (1982).
  * Bernard [1983] C. W. Bernard, Nucl. Phys. B 219, 341 (1983).
  * Donoghue [1984] J. F. Donoghue, Phys. Rev. D 29, 2559 (1984).
  * Mandula and Ogilvie [1987] J. Mandula and M. Ogilvie, Phys. Lett. B 185, 127 (1987).
  * Cornwall and Papavassiliou [1989] J. M. Cornwall and J. Papavassiliou, Phys. Rev. D 40, 3474 (1989).
  * Lavelle [1991] M. Lavelle, Phys. Rev. D 44, 26 (1991).
  * Halzen  _et al._ [1993] F. Halzen, G. I. Krein, and A. A. Natale, Phys. Rev. D 47, 295 (1993).
  * 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).
  * Mihara and Natale [2000] A. Mihara and A. A. Natale, Phys. Lett. B 482, 378 (2000).
  * Kondo [2001] K.-I. Kondo, Phys. Lett. B 514, 335 (2001).
  * Philipsen [2002] O. Philipsen, Nucl. Phys. B628, 167 (2002).
  * Aguilar  _et al._ [2003] A. C. Aguilar, A. A. Natale, and P. S. Rodrigues da Silva, Phys. Rev. Lett. 90, 152001 (2003).
  * Pawlowski  _et al._ [2004] J. M. Pawlowski, D. F. Litim, S. Nedelko, and L. von Smekal, Phys. Rev. Lett. 93, 152002 (2004).
  * Aguilar and Natale [2004] A. C. Aguilar and A. A. Natale, JHEP 08, 057 (2004).
  * Aguilar and Papavassiliou [2006] A. C. Aguilar and J. Papavassiliou, JHEP 12, 012 (2006).
  * Epple  _et al._ [2008] D. Epple, H. Reinhardt, W. Schleifenbaum, and A. P. Szczepaniak, Phys. Rev. D77, 085007 (2008).
  * Aguilar  _et al._ [2008] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D78, 025010 (2008).
  * Braun  _et al._ [2010] J. Braun, H. Gies, and J. M. Pawlowski, Phys. Lett. B684, 262 (2010).
  * Tissier and Wschebor [2010] M. Tissier and N. Wschebor, Phys. Rev. D 82, 101701 (2010).
  * Campagnari and Reinhardt [2010] D. R. Campagnari and H. Reinhardt, Phys. Rev. D82, 105021 (2010).
  * Tissier and Wschebor [2011] M. Tissier and N. Wschebor, Phys. Rev. D84, 045018 (2011).
  * Serreau and Tissier [2012] J. Serreau and M. Tissier, Phys. Lett. B712, 97 (2012).
  * Fagundes  _et al._ [2012] D. A. Fagundes, E. G. S. Luna, M. J. Menon, and A. A. Natale, Nucl. Phys. A 886, 48 (2012).
  * Binosi  _et al._ [2012] D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D86, 085033 (2012).
  * Maas [2013] A. Maas, Phys. Rept. 524, 203 (2013).
  * 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).
  * Cyrol  _et al._ [2018] A. K. Cyrol, J. M. Pawlowski, A. Rothkopf, and N. Wink, SciPost Phys. 5, 065 (2018).
  * Gracey  _et al._ [2019a] J. A. Gracey, M. Peláez, U. Reinosa, and M. Tissier, Phys. Rev. D 100, 034023 (2019a).
  * Roberts [2020] C. D. Roberts, Symmetry 12, 1468 (2020).
  * Peláez  _et al._ [2021] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, Rept. Prog. Phys. 84, 124202 (2021).
  * Eichmann  _et al._ [2021] G. Eichmann, J. M. Pawlowski, and J. a. M. Silva, Phys. Rev. D 104, 114016 (2021).
  * Horak  _et al._ [2022] J. Horak, F. Ihssen, J. Papavassiliou, J. M. Pawlowski, A. Weber, and C. Wetterich, SciPost Phys. 13, 042 (2022).
  * Papavassiliou [2022] J. Papavassiliou, Chin. Phys. C 46, 112001 (2022).
  * Ding  _et al._ [2023] M. Ding, C. D. Roberts, and S. M. Schmidt, Particles 6, 57 (2023).
  * Ferreira and Papavassiliou [2023] M. N. Ferreira and J. Papavassiliou, Particles 6, 312 (2023).
  * Poggio  _et al._ [1975] E. C. Poggio, E. Tomboulis, and S. H. H. Tye, Phys. Rev. D11, 2839 (1975).
  * Aguilar  _et al._ [2012] A. C. Aguilar, D. Ibanez, V. Mathieu, and J. Papavassiliou, Phys. Rev. D85, 014018 (2012).
  * Ibañez and Papavassiliou [2013] D. Ibañez and J. Papavassiliou, Phys. Rev. D87, 034008 (2013).
  * Aguilar  _et al._ [2016a] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Front. Phys.(Beijing) 11, 111203 (2016a).
  * Aguilar  _et al._ [2018] A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Eur. Phys. J. C78, 181 (2018).
  * Aguilar  _et al._ [2022] A. C. Aguilar, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 105, 014030 (2022).
  * Roberts and Williams [1994] C. D. Roberts and A. G. Williams, Prog. Part. Nucl. Phys. 33, 477 (1994).
  * Alkofer and von Smekal [2001] R. Alkofer and L. von Smekal, Phys. Rept. 353, 281 (2001).
  * Maris and Roberts [2003] P. Maris and C. D. Roberts, Int. J. Mod. Phys. E12, 297 (2003).
  * Fischer [2006] C. S. Fischer, J. Phys. G 32, R253 (2006).
  * Roberts [2008] C. D. Roberts, Prog. Part. Nucl. Phys. 61, 50 (2008).
  * Fischer  _et al._ [2009] C. S. Fischer, A. Maas, and J. M. Pawlowski, Annals Phys. 324, 2408 (2009).
  * Binosi and Papavassiliou [2009] D. Binosi and J. Papavassiliou, Phys. Rept. 479, 1 (2009).
  * 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).
  * Fister and Pawlowski [2013] L. Fister and J. M. Pawlowski, Phys. Rev. D88, 045010 (2013).
  * Cloet and Roberts [2014] I. C. Cloet and C. D. Roberts, Prog. Part. Nucl. Phys. 77, 1 (2014).
  * Binosi  _et al._ [2015] D. Binosi, L. Chang, J. Papavassiliou, and C. D. Roberts, Phys. Lett. B742, 183 (2015).
  * Binosi  _et al._ [2016] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D93, 096010 (2016).
  * Binosi  _et al._ [2017] D. Binosi, C. Mezrag, J. Papavassiliou, C. D. Roberts, and J. Rodriguez-Quintero, Phys. Rev. D96, 054026 (2017).
  * Huber [2020a] M. Q. Huber, Phys. Rept. 879, 1 (2020a).
  * Huber [2020b] M. Q. Huber, Phys. Rev. D 101, 114009 (2020b).
  * Aguilar  _et al._ [2016b] A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D94, 045002 (2016b).
  * 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).
  * Alexandrou  _et al._ [2002] C. Alexandrou, P. de Forcrand, and E. Follana, Phys. Rev. D 65, 114508 (2002).
  * Cucchieri and Mendes [2007] A. Cucchieri and T. Mendes, PoS LATTICE2007, 297 (2007).
  * Bogolubsky  _et al._ [2007] I. Bogolubsky, E. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, PoS LATTICE2007, 290 (2007).
  * Bowman  _et al._ [2007] P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, A. Sternbeck, L. von Smekal, A. G. Williams, and J.-b. Zhang, Phys. Rev. D 76, 094505 (2007).
  * Kamleh  _et al._ [2007] W. Kamleh, P. O. Bowman, D. B. Leinweber, A. G. Williams, and J. Zhang, Phys. Rev. D76, 094501 (2007).
  * Cucchieri and Mendes [2008] A. Cucchieri and T. Mendes, Phys. Rev. Lett. 100, 241601 (2008).
  * Bogolubsky  _et al._ [2009] I. Bogolubsky, E. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Phys. Lett. B676, 69 (2009).
  * Oliveira and Silva [2009] O. Oliveira and P. Silva, PoS LAT2009, 226 (2009).
  * Cucchieri  _et al._ [2009] A. Cucchieri, T. Mendes, and E. M. S. Santos, Phys. Rev. Lett. 103, 141602 (2009).
  * Cucchieri and Mendes [2010] A. Cucchieri and T. Mendes, Phys. Rev. D81, 016005 (2010).
  * Cucchieri  _et al._ [2010] A. Cucchieri, T. Mendes, G. M. Nakamura, and E. M. S. Santos, PoS FACESQCD, 026 (2010).
  * Oliveira and Bicudo [2011] O. Oliveira and P. Bicudo, J. Phys. G G38, 045003 (2011).
  * Ayala  _et al._ [2012] A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, and J. Rodriguez-Quintero, Phys. Rev. D86, 074512 (2012).
  * Sternbeck and Müller-Preussker [2013] A. Sternbeck and M. Müller-Preussker, Phys. Lett. B 726, 396 (2013).
  * Bicudo  _et al._ [2015] P. Bicudo, D. Binosi, N. Cardoso, O. Oliveira, and P. J. Silva, Phys. Rev. D92, 114514 (2015).
  * Duarte  _et al._ [2016] A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 014502 (2016).
  * Dudal  _et al._ [2018] D. Dudal, O. Oliveira, and P. J. Silva, Annals Phys. 397, 351 (2018).
  * 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._ [2023b] A. C. Aguilar, M. . N. Ferreira, D. Ibañez, and J. Papavassiliou, Eur. Phys. J. C 83, 967 (2023b).
  * Vladimirov [1971] V. Vladimirov,  _Equations of Mathematical Physics_, Monographs and textbooks in pure and applied mathematics (M. Dekker, 1971).
  * Polyanin and Manzhirov [2008] A. Polyanin and A. Manzhirov,  _Handbook of Integral Equations: Second Edition_ (Taylor & Francis, 2008).
  * 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,  _Proceedings, Laws Of Hadronic Matter, Erice,_ (MIT, Cambridge, MA, 1973).
  * Cornwall and Norton [1973] J. Cornwall and R. Norton, Phys. Rev. D 8, 3338 (1973).
  * 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).
  * Boucaud  _et al._ [2017] P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D95, 114503 (2017).
  * 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).
  * 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).
  * Maas and Vujinović [2022] A. Maas and M. Vujinović, SciPost Phys. Core 5, 019 (2022).
  * Catumba  _et al._ [2022a] G. T. R. Catumba, O. Oliveira, and P. J. Silva, PoS LATTICE2021, 467 (2022a).
  * Catumba  _et al._ [2022b] G. T. R. Catumba, O. Oliveira, and P. J. Silva, EPJ Web Conf. 258, 02008 (2022b).
  * 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).
  * Berges [2004] J. Berges, Phys. Rev. D 70, 105010 (2004).
  * Carrington and Guo [2011] M. E. Carrington and Y. Guo, Phys. Rev. D 83, 016006 (2011).
  * York  _et al._ [2012] M. C. A. York, G. D. Moore, and M. Tassler, JHEP 06, 077 (2012).
  * Mueller and Pawlowski [2015] N. Mueller and J. M. Pawlowski, Phys. Rev. D91, 116010 (2015).
  * Williams  _et al._ [2016] R. Williams, C. S. Fischer, and W. Heupel, Phys. Rev. D93, 034026 (2016).
  * Aguilar  _et al._ [2023c] A. C. Aguilar, M. N. Ferreira, J. Papavassiliou, and L. R. Santos, Eur. Phys. J. C 83, 549 (2023c).
  * Huber  _et al._ [2020] M. Q. Huber, C. S. Fischer, and H. Sanchis-Alepuz, Eur. Phys. J. C 80, 1077 (2020).
  * Huber  _et al._ [2021] M. Q. Huber, C. S. Fischer, and H. Sanchis-Alepuz, Eur. Phys. J. C 81, 1083 (2021), [Erratum: Eur.Phys.J.C 82, 38 (2022)].
  * Taylor [1971] J. Taylor, Nucl. Phys. B 33, 436 (1971).
  * Slavnov [1972] A. Slavnov, Theor. Math. Phys. 10, 99 (1972).
  * Aguilar  _et al._ [2023d] A. C. Aguilar, M. N. Ferreira, B. M. Oliveira, J. Papavassiliou, and L. R. Santos, Eur. Phys. J. C 83, 889 (2023d).
  * Ball and Chiu [1980] J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2550 (1980), [Erratum: Phys. Rev. D 23, 3085 (1981)].
  * 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)].
  * Gracey  _et al._ [2019b] J. Gracey, H. Kißler, and D. Kreimer, Phys. Rev. D 100, 085001 (2019b).
  * Aguilar  _et al._ [2014] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D89, 085032 (2014).
  * Bjorken and Drell [1965] J. D. Bjorken and S. D. Drell,  _Relativistic quantum fields_ , International Series In Pure and Applied Physics (McGraw-Hill, New York, 1965).
  * Cui  _et al._ [2020] Z.-F. Cui, J.-L. Zhang, D. Binosi, F. de Soto, C. Mezrag, J. Papavassiliou, C. D. Roberts, J. Rodríguez-Quintero, J. Segovia, and S. Zafeiropoulos, Chin. Phys. C 44, 083102 (2020).
  * 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).
  * Mitter  _et al._ [2015] M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D91, 054035 (2015).
  * Pawlowski  _et al._ [2022] J. M. Pawlowski, C. S. Schneider, and N. Wink, arXiv:2202.11123 [hep-th] (2022).
  * Press  _et al._ [1992] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in FORTRAN: The Art of Scientific Computing_ (Cambridge University Press, 1992).
  * Altarelli [1982] G. Altarelli, Phys. Rept. 81, 1 (1982).
  * von Smekal  _et al._ [1998] L. von Smekal, A. Hauck, and R. Alkofer, Annals Phys. 267, 1 (1998), [Erratum: Annals Phys. 269, 182 (1998)].
  * Fischer  _et al._ [2002] C. S. Fischer, R. Alkofer, and H. Reinhardt, Phys. Rev. D 65, 094008 (2002).
  * Pennington and Wilson [2011] M. R. Pennington and D. J. Wilson, Phys. Rev. D84, 119901 (2011).
  * Celmaster and Gonsalves [1979a] W. Celmaster and R. J. Gonsalves, Phys. Rev. Lett. 42, 1435 (1979a).
  * Celmaster and Gonsalves [1979b] W. Celmaster and R. J. Gonsalves, Phys. Rev. D 20, 1420 (1979b).
  * Alles  _et al._ [1997] B. Alles, D. Henty, H. Panagopoulos, C. Parrinello, C. Pittori, and D. G. Richards, Nucl. Phys. B502, 325 (1997).
  * Boucaud  _et al._ [2009] P. Boucaud, F. De Soto, J. Leroy, A. Le Yaouanc, J. Micheli, _et al._ , Phys. Rev. D79, 014508 (2009).
  * Huber [2017] M. Q. Huber, Eur. Phys. J. C77, 733 (2017).
  * Binosi  _et al._ [2014] D. Binosi, D. Ibañez, and J. Papavassiliou, JHEP 09, 059 (2014).
  * Cyrol  _et al._ [2015] A. K. Cyrol, M. Q. Huber, and L. von Smekal, Eur. Phys. J. C75, 102 (2015).
  * Pawlowski  _et al._ [2023] J. M. Pawlowski, C. S. Schneider, J. Turnwald, J. M. Urban, and N. Wink, Phys. Rev. D 108, 076018 (2023).
  * Colaço  _et al._ [2024] M. Colaço, O. Oliveira, and P. J. Silva, Phys. Rev. D 109, 074502 (2024).
  * Aguilar  _et al._ [2024] A. C. Aguilar, M. N. Ferreira, J. Papavassiliou, and L. R. Santos, arXiv:2402.16071 [hep-ph] (2024).