The soft-gluon limit and the infrared enhancement of the quark-gluon vertex

Abstract

The Schwinger–Dyson quark equation (SDE) combined with results from lattice simulation for the propagators are used to obtain information on the quark-gluon vertex, taking into account the recent full QCD lattice results for the soft-gluon limit. Its inclusion leads to a clear enhancement of the infrared quark-gluon vertex. We also find that the relative contribution of the quark-ghost kernel to the quark-gluon vertex in the infrared region does not follow the rules from the perturbative analysis of the ultraviolet region. This shows that for QCD the intuition based on perturbation theory does not apply to the full momentum range. The framework developed in the current work provides analytical expressions for all the longitudinal components of vertex taken into account.

1 Introduction and motivation

The interaction between quarks and gluons is described by Quantum Chromodynamics (QCD). If at high energies asymptotic freedom allows the use of perturbation theory, at low energies QCD becomes strongly coupled and other type of techniques are needed to solve the theory. Confinement and chiral symmetry breaking are two important non-perturbative outcomes of QCD associated with its infrared (IR) properties [1, 2]. In particular, confinement impacts on the computation of hadronic transitions elements also at high energies.

Lattice QCD simulations is a successful ab initio approach to study the non-perturbative regime. The use of continuous methods, such as Schwinger–Dyson equations (SDE), Bethe–Salpeter, Faddeev equations, functional renormalisation group techniques and Holographic models complement our catalogue of non-perturbative tools for quantum field theories.

The gluon, the ghost and the quark propagators have been thoroughly studied within the non-perturbative regime of QCD and, nowadays, one has a fair understanding and description of these two point correlation functions. On the other hand, our level of knowledge on the three particle vertices is much poor when compared to the two point functions. Among the QCD three particle vertices is the quark-gluon vertex that plays a major role in our understanding of hadrons.

Herein we pursuit the efforts initiated in [3] and developed further in [4, 5] that combine results from lattice simulation in the Landau gauge, for the propagators, and the quark gap equation, together with an exact solution of the Slavnov–Taylor identity for the vertex, to get information on the quark-gluon vertex. The exact solution of the Slavnov-Taylor identity for the vertex calls for the quark-ghost kernel, a four point Green function, that is described in terms of four form factors [6], our unknowns. The main novelty in the present work is the inclusion of the lattice results for the soft-gluon limit of the quark-gluon vertex, that translate into normalisation conditions for the quark-ghost kernel form factors. In practice, this was achieved by building an ansatz that is inspired on the analysis of the lattice data for the quark-gluon vertex [4] and the gap equation is solved for the form factors, written as Padé approximants. We stress that all other quantities in the gap equation are known from lattice simulations, including the quark propagator in the Landau gauge. The resulting quark-gluon form factors are identical to those computed previously [5] with the exception of the infrared region, where it is observed a clear enhancement. Furthermore, a detailed analysis of the relative contribution of the quark-ghost kernel form factors to the quark-gluon vertex does not comply, in the infrared region, with their ultraviolet (perturbative) behaviour.

This paper reports on results for the quark-gluon vertex as outlined above and it is organized as follows. In Sect. 2, we set our notation for the quark gap equation and the quark-gluon vertex, where as usual we consider only the longitudinal components, written in terms of the quark-ghost kernel form factors, which are parametrized by Padé approximants constrained by the soft-gluon limit results from lattice simulations. In Sect. 3, the Schwinger–Dyson equation is decomposed for the vector and scalar parts of the quark self energies, and the minimal requirements for the ultraviolet behaviour of the quark-ghost kernel are discussed together with its anzats. Also, the explicit form of our anzatsë for the longitudinal form factors are provided. In Sect. 4, we built our anzats such that the soft-gluon limit from QCD lattice simulations of $\lambda _1$ is incorporated. In Sect. 5, we present the results from the inversion of the Schwinger–Dyson equations to get the coefficients of the Padé approximants for the quark-ghost kernel relying on simulating annealing to minimize the sum of the relative error of the scalar and vector equations. In Sect. 6 we present our results for the form factors of the longitudinal components of the quark-gluon vertex and also analyze their contribution of the quark-ghost kernel separately. In Sect. 7 we provide a discussion of our results to put in the perspective of some previous studies. This section is closed with a summary of our work.

2 The Quark gap equation and the Quark-Gluon Vertex

The quark propagator is color diagonal and its spin-Lorentz structure reads, in Minkowski space,

$$\begin{aligned} S^{-1} (p)= & {} -i \, \left( A(p^2) \,/\!\! p\ - B(p^2) \right) \nonumber \\= & {} -i \, Z_2 (\, \,/\!\! p - m^{\mathrm {bm}}) + \varSigma (p) \end{aligned}$$

(1)

where $Z(p^2) = 1/A(p^2)$ is the quark wave function, $M(p^2) = B(p^2) / A(p^2)$ the renormalisation group invariant running quark mass, $Z_2$ is the quark renormalisation constant and $m^{\mathrm {bm}}$ the bare current quark mass. The quark self-energy is given by

$$\begin{aligned} \varSigma (p)= & {} Z_1 \, \int {d^4 q\over (2\pi )^4} ~ D_{\mu \nu }^{a b}(q) ~ (\, i\, g\, t^b \gamma _\nu \,) ~ \nonumber \\&\qquad \quad \times S(p-q) ~ \varGamma ^a_\mu (-p, \, p-q, \, q) , \end{aligned}$$

(2)

where $Z_1$ is a combination of several renormalisation constants and the Landau gauge gluon propagator is

$$\begin{aligned} D^{ab}_{\mu \nu } (q) = -i\, \delta ^{ab} \left( g_{\mu \nu } - \frac{q_\mu q_\nu }{q^2} \right) D ( q^2 ) \ . \end{aligned}$$

(3)

The quark-gluon vertex is defined with incoming momenta $p_1+p_2+p_3 = 0$, where $p_2$ is the incoming quark momentum, $-p_1$ the outgoing quark momentum and $p_3$ the incoming gluon momentum. Our notation follows that used in [4, 6]. The one-particle irreducible quark-gluon Green function is depicted as

$$\begin{aligned} \varGamma ^a_\mu (p_1, p_2, p_3) = g \, t^a \, \varGamma _\mu (p_1, p_2, p_3) \ , \end{aligned}$$

(4)

where g is the strong coupling constant and $t^a$ are the generators of the color SU(3) group in the fundamental representation.

Assuming that the gluon propagator and $\varGamma ^a_\mu$ are known, from the gap equation (2), one can get the quark propagator. If $Z(p^2)$ and $M(p^2)$ are known, it is possible to use (2) to extract information on the quark-gluon vertex. From the mathematical point of view, computing $\varGamma ^a_\mu$ from the gap equation means solving an ill-defined problem. The introduction of a prior, that can be accommodated by regularising the integral equation or introducing a basis of functions, allows to exactly and unambiguously solve the modified equation for the vertex. The solution depends on the prior and one should check its (in)dependence on the prior.

The vertex function $\varGamma _\mu$, see Eq. (4), can be decomposed in a longitudinal $\varGamma ^{(L)}_\mu$ and a transverse $\varGamma ^{(T)}_\mu$ component, relative to the gluon momenta, as

$$\begin{aligned} \varGamma _{\mu }(p_1, \, p_2, \, p_3) = \varGamma ^{(L)}_{\mu }(p_1, \, p_2, \, p_3) + \varGamma ^{(T)}_{\mu }(p_1, \, p_2, \, p_3) \end{aligned}$$

(5)

and, by definition,

$$\begin{aligned} p_{3}^{\mu } ~ \varGamma ^{(T)}_{\mu }(p_1, \, p_2, \, p_3)= 0 \ . \end{aligned}$$

(6)

As is usual in the analysis of the Dyson-Schwinger equations, in the current work we will focus on the longitudinal component of the $\varGamma _\mu$ and will ignore $\varGamma ^{(T)}_{\mu }$. If a tensor basis for $\varGamma ^{(L)}_\mu$ and $\varGamma ^{(T)}_\mu$ is given, then $\varGamma _\mu$ is a sum of scalar form factors that multiply each of the elements of the tensor basis. The full vertex requires twelve form factors, with four of them being associated with the longitudinal component and eight with the transverse one. In the Ball-Chiu (BC) basis [7], the longitudinal vertex is written as:

$$\begin{aligned} \varGamma ^\mathrm {L}_\mu (p_1, p_2, p_3)= & {} -i \, \bigg ( \lambda _1 \, \gamma _\mu ~ + ~ \lambda _2 \, ( \,/\!\! p_1 - \,/\!\! p_2 ) \left( p_1 - p_2 \right) _\mu \nonumber \\&+ ~ \lambda _3 \, ( p_1 - p_2 )_\mu ~ + ~ \lambda _4 \, \sigma _{\mu \nu } \left( p_1 - p_2 \right) ^\nu \bigg ) , \nonumber \\ \, \end{aligned}$$

(7)

where $\sigma _{\mu \nu } = {1\over 2} [\gamma _{\mu },\gamma _{\nu }]$ and $\lambda _i = \lambda _i(p^2_1, p^2_2,p^2_3)$. The transverse component of the BC vertex is written as

$$\begin{aligned} \varGamma ^\mathrm {T}_\mu (p_1, p_2, p_3)= & {} -i\, \sum ^8_{i=1} \tau _i (p_1, p_2, p_3) ~ T^{(i)}_\mu (p_1 , p_2) \ , \end{aligned}$$

(8)

where the operators associated to the transverse part of the vertex are

$$\begin{aligned} T^{(1)}_\mu (p_1 , p_2)= & {} \left[ p_{1\mu } \left( p_2 \cdot p_3 \right) - p_{2\mu } \left( p_1 \cdot p_3 \right) \right] {\mathbb {I}}_D \ , \nonumber \\ T^{(2)}_\mu (p_1 , p_2)= & {} - T^{(1)}_\mu (p_1 , p_2)~ \left( \,/\!\! p_1 - \,/\!\! p_2 \right) \ , \nonumber \\ T^{(3)}_\mu (p_1 , p_2)= & {} p^2_3 \, \gamma _\mu - p_{3\mu } \ \,/\!\! p_3 \ , \nonumber \\ T^{(4)}_\mu (p_1 , p_2)= & {} T^{(1)}_\mu (p_1 , p_2) ~ \sigma _{\alpha \beta }\, p^\alpha _1 p^\beta _2 \ , \nonumber \\ T^{(5)}_\mu (p_1 , p_2)= & {} \sigma _{\mu \nu } \, p^\nu _3 \ , \nonumber \\ T^{(6)}_\mu (p_1 , p_2)= & {} \gamma _\mu \left( p^2_1 - p^2_2 \right) + \left( p_1 - p_2 \right) _\mu \, \,/\!\! p_3 \ , \nonumber \\ T^{(7)}_\mu (p_1 , p_2)= & {} -{1\over 2} \left( p^2_1 - p^2_2 \right) \left[ \gamma _\mu \left( \,/\!\! p_1 - \,/\!\! p_2 \right) - \left( p_1 - p_2 \right) _\mu {\mathbb {I}}_D \right] \nonumber \\&\quad - \left( p_1 - p_2 \right) _\mu \sigma _{\alpha \beta } \, p^\alpha _1 \, p^\beta _2 \ , \nonumber \\ T^{(8)}_\mu (p_1 , p_2)= & {} - \gamma _\mu \, \sigma _{\alpha \beta } \, p^\alpha _1 \, p^\beta _2 + ( p_{1\mu } \, \,/\!\! p_2 - p_{2\mu } \, \,/\!\! p_1) \ . \end{aligned}$$

(9)

The discrete symmetries of QCD constraint the quark-gluon vertex. For the longitudinal form factors, charge conjugation invariance [6] requires

$$\begin{aligned} \lambda _i \left( p^2_1, \, p^2_2, \, p^2_3\right)= & {} \lambda _i \left( p^2_2, \, p^2_1, \, p^2_3\right) |_{i=1,2,3}, \nonumber \\ \lambda _4 \left( p^2_1, \, p^2_2, \, p^2_3\right)= & {} - \lambda _4 \left( p^2_2, \, p^2_1, \, p^2_3\right) \ . \end{aligned}$$

(10)

These properties under interchange of quark momenta imply that when $p^2_1 = p^2_2$ the form factor $\lambda _4 = 0$ as happens for the soft-gluon limit where $p_3 = 0$.

The Slavnov-Taylor identity (STI) for the quark-gluon vertex,

$$\begin{aligned} p^\mu _3 ~ \varGamma _\mu (p_1, \, p_2, \, p_3 )= & {} F(p^2_3) \Big [ S^{-1}(-p_1) \, H(p_1, \, p_2, \, p_3)\nonumber \\&- {\overline{H}}(p_2, \, p_1, \, p_3) \, S^{-1}(p_2) \Big ] , \end{aligned}$$

(11)

relates the longitudinal vertex form factors, $\lambda _i$’s $(i=1-4)$, with the quark propagator, the quark-ghost kernels that are define by H and ${\overline{H}}$, see [6] for notation and definitions, and the ghost dressing function $F(p^2_3)$, related to the ghost propagator by

$$\begin{aligned} D^{ab}(p^2) = - i \, \delta ^{ab} \, F(p^2) / p^2 \ . \end{aligned}$$

(12)

The quark-ghost kernel can be written in terms of four form factors [6], called $X_i\equiv X_i(p^2_1, p^2_2, p^2_3)$ and ${\overline{X}}_i \equiv X_i (p^2_2, \, p^2_1, \, p^2_3 )$ , as

$$\begin{aligned} H(p_1, \, p_2, \, p_3 )= & {} X_0 {\mathbb {I}}_D + X_1 \, / \!\!p_1 + X_2 \, / \!\!p_2 + X_3 \,\sigma _{\alpha \beta } \, p^\alpha _1 \, p^\beta _2 \ , \nonumber \\ {\overline{H}}(p_2, \, p_1, \, p_3 )= & {} X_0 {\mathbb {I}}_D - X_2 \, / \!\!p_1 - X_1 \, /\!\!p_2 + X_3 \,\sigma _{\alpha \beta } \, p^\alpha _1 \, p^\beta _2 \ .\nonumber \\ \end{aligned}$$

(13)

The STI (11) can be solved exactly for the $\lambda _i$’s [8] resulting in

$$\begin{aligned} \lambda _1= & {} \frac{F(p^2_3)}{2} \Bigg \{ A(p^2_1) \left[ X_0 + \left( p^2_1 - p_1 \cdot p_2 \right) X_3 \right] \nonumber \\&+ B(p^2_1) \left[ X_1 + X_2 \right] \Bigg \} + (p_1\leftrightarrow p_2) \ , \nonumber \\ \lambda _2= & {} \frac{F(p^2_3)}{2 (p^2_2 - p^2_1)} \Bigg \{ A(p^2_1)\! \left[ \left( p^2_1 + p_1\! \cdot p_2 \right) X_3 - X_0 \right] \nonumber \\&+ B(p^2_1) \left[ X_2 - X_1 \right] \Bigg \}+ (p_1\leftrightarrow p_2) \ , \nonumber \\ \lambda _3= & {} \frac{F(p^2_3)}{p^2_1 - p^2_2} \Bigg \{ A(p^2_1) \left[ p^2_1 \, X_1 + p_1 \cdot p_2 \ X_2\right] + B(p^2_1) \, X_0 \Bigg \} \nonumber \\ \nonumber \\ \lambda _4= & {} - \frac{F(p^2_3)}{2} \Bigg \{ A(p^2_1) \, X_2 + B(p^2_1) \, X_3 \Bigg \} - (p_1\leftrightarrow p_2) \ . \nonumber \\ \end{aligned}$$

(14)

These expressions comply with the symmetries under interchange of quark momenta listed in (10). Moreover, as discussed in [5] all $\lambda _i$ are regular functions.

The tree level perturbative solution of QCD for the quark-ghost kernel form factors gives $X_0 = 1$ and $X_{1,2,3} = 0$ [6] and, in this case, the longitudinal quark-gluon vertex is uniquely determined by the quark wave function $Z(p^2)$ and the running quark mass $M(p^2)$. On the other hand, the non-perturbative solutions of QCD [4, 9,10,11,12] result in $X_{1,2,3}$ that deviate significantly from their tree level values, specially at infrared mass scales. These non-perturbative solutions for $X_0$ return a function that hardly deviates from its tree level perturbative value [3, 4, 9, 12, 13] and $X_0 \approx 1$ seems to be a good approximation for this form factor. We call the reader’s attention that so far we only have approximate non-perturbative solutions for QCD. At a qualitative level, the various approaches seem to produce compatible results. See also [14,15,16,17] for other recent studies of the quark-gluon vertex.

There are now results from full QCD lattice simulations for $\lambda _1(p^2)$ in the soft-gluon limit [18]. This information should be incorporated in the inversion of the gap equation to compute the quark-gluon vertex. The longitudinal form factors $\lambda _1$ in this limit [5] is given by

$$\begin{aligned} \lambda _1 (p^2) = \frac{F(0)}{Z(p^2)} \bigg \{ 1 ~ + ~ 2 \, X_1(p^2) \, M(p^2) ~ + ~ 2 \, p^2 \, X_3 (p^2) \bigg \} \ , \end{aligned}$$

(15)

where p is the quark momentum. This kinematical configuration is interesting also because it gets no contribution from the transverse form factors. This can be easily checked using the transverse basis considered in [19].

3 Decomposing the Schwinger–Dyson equation

In the current study we use only the longitudinal components of the quark-gluon vertex; see Eq. (7) and Sec. 3.2 for detailed expressions for all the form factors. By taking traces of (1) one accesses either its scalar or vector part. In order to include inputs from lattice simulations, we perform a Wick rotation in the Schwinger–Dyson equations, writing them in the Euclidean space (see e.g. [4] and references therein) as

$$\begin{aligned}&B(p^2) = Z_2\, m^\mathrm {bm} \nonumber \\&\qquad + ~ C_F Z_1 g^2 \, \int \frac{d^4 q}{(2 \, \pi )^4} ~ \frac{ Z(k^2) \, D(q^2) }{ k^2 + \left[ M(k^2) \right] ^2} \Bigg \{ {\mathcal {A}}_s + M(k^2) \, {\mathcal {B}}_s \Bigg \} , \end{aligned}$$

(16)

$$\begin{aligned}&p^2 A(p^2) = Z_2\, p^2 \nonumber \\&\qquad - ~ C_F Z_1 g^2 \int \frac{d^4 q}{(2 \, \pi )^4} ~ \frac{ Z(k^2) \, D(q^2) }{ k^2 + \left[ M(k^2) \right] ^2} \Bigg \{ {\mathcal {A}}_v + M(k^2) \, {\mathcal {B}}_v \, \Bigg \}\nonumber \\ \end{aligned}$$

(17)

with

$$\begin{aligned}&{\mathcal {A}}_s(p^2,k^2,q^2) = 2 ~\lambda _3 ~ \left( p^2 - \frac{(pq)^2}{q^2} \right) \nonumber \\&\quad + ~ \lambda _4 ~ \left[ 9 \, (pq) - 2 \, \frac{(pq)^2}{q^2} - 4 \, p^2 - 3 \, q^2 \right] \, , \end{aligned}$$

(18)

$$\begin{aligned}&{\mathcal {B}}_s(p^2,k^2,q^2) = 3 ~\lambda _1 + 4 ~\lambda _2 ~\left[ \frac{(pq)^2}{q^2} - p^2 \right] \, , \end{aligned}$$

(19)

$$\begin{aligned}&{\mathcal {A}}_v(p^2,k^2,q^2) = \lambda _1 ~ \left[ p^2 + 2 \, \frac{(pq)^2}{q^2} - 3 \, (pq) \right] \nonumber \\&\quad + ~2 \, \lambda _2 ~ \bigg [ p^2 q^2 + 2 \, \frac{(pq)^3}{q^2} \nonumber \\&\quad - ~\left( 2 \, \frac{p^2}{q^2} + 1 \right) (pq)^2 - 2 \, p^2 (pq) + 2 \, p^4 \bigg ] , \end{aligned}$$

(20)

$$\begin{aligned}&{\mathcal {B}}_v(p^2,k^2,q^2) = 2 ~ \lambda _3 \left[ p^2 - \frac{(pq)^2}{q^2} \right] \nonumber \\&\quad + ~ \lambda _4 \left[ 4 \, p^2 + 2 \, \frac{(pq)^2}{q^2} - 3 \, (pq) \right] , \end{aligned}$$

(21)

where $k = p -q$, $C_F = 4/3$ is the Casimir invariant associated to the SU(3) fundamental representation and we used $\lambda _i\equiv \lambda _i (p^2 , \, k^2, \, q^2)$ to simplify the notation.

3.1 Results from 1-loop perturbation theory

The full set of $\lambda _i$’s were computed in [20] in perturbation theory to one-loop for arbitrary linear covariant gauges. The authors report full expressions for the various form factors at several kinematic configurations. For the symmetric limit where $p^2 = k^2 = q^2$, the leading behaviour for the non-vanishing form factors are

$$\begin{aligned} \lambda _1 (p^2) ~ \propto ~ 1 \ , \qquad \lambda _2 (p^2)\,\, \text {and} \,\,\lambda _3 (p^2) ~ \propto ~ 1 / p^2 \end{aligned}$$

(22)

up to logarithmic corrections. The solution of the STI for $\lambda _1$ reads

$$\begin{aligned}&\lambda _1 (p^2) = F(p^2) \Bigg \{ A(p^2) \bigg [ X_0 + \frac{5}{2} \, p^2 \, X_3 \bigg ] + 2 \, B(p^2) \, X_1 \Bigg \} \ \nonumber \\&\quad \hbox { and }~ X_2 = X_1 \, . \end{aligned}$$

(23)

At large quark momentum $\lambda _1 (p^2)$, $F(p^2)$, $A(p^2)$ and $B(p^2)$ become constants and this expression suggests that $X_0 \sim 1$, $X_1 \sim 1$ and $X_3 \sim 1/p^2$. These are minimal scaling laws in the sense that some of the form factors can behave as larger negative power of $p^2$. A similar analysis of the symmetric limit for $\lambda _2$ and $\lambda _3$ returns minimal scaling laws compatible with the above ones. The tree level expressions for the quark-ghost kernel form factors are $X_0 = 1$ and $X_{1, 2, 3} = 0$ and these results should be recovered in the UV limit. This suggests that the minimal scaling law for $X_1$ is not that previously reported but, instead, a larger negative power of $p^2$ as e.g.

$$\begin{aligned} X_0 \sim 1 ~ , \qquad X_1 \sim 1/p^2 \qquad \hbox { and }\qquad X_3 \sim 1/p^2 \ . \end{aligned}$$

(24)

This results are in good agreement with the scaling analysis performed in [4] and are compatible with the minimal scaling laws derived directly from one-loop perturbation theory quoted above.

3.2 The longitudinal Quark-Gluon vertex

In order to compute a solution for vertex from the SDE we assume functional forms for the $X_i$’s. As in [4] the gluon, the ghost [21] and the quark [22] propagators are taken from lattice simulations and modelled to reproduce its correct UV perturbative behaviour. Exact expressions can be found in App. A. Moreover, as discussed in [4], following the analysis of the lattice data for soft-gluon limit performed in [5], an ansatz for the quark-ghost form factors is

$$\begin{aligned} X_0 (p^2, k^2, q^2)\equiv & {} X_0(q^2) , \end{aligned}$$

(25)

$$\begin{aligned} X_1 (p^2, k^2, q^2)= & {} X_2 (p^2, k^2, q^2) \equiv D\left( \frac{p^2 + k^2}{2} \right) ~ Y_1(q^2) , \nonumber \\ \end{aligned}$$

(26)

$$\begin{aligned} X_3 (p^2, k^2, q^2)\equiv & {} D\left( \frac{p^2 + k^2}{2} \right) ~ Y_3(q^2) \ . \end{aligned}$$

(27)

This ansatz solves the SDE with a relative error smaller than 4% - see [4] for further details.

The parameterisations (25) - (27) have to comply with the exact results derived in [9]

$$\begin{aligned} X_0 (p^2, \, p^2, \, 0) = 1 ~ \hbox { and } ~ X_1(p^2, \, p^2, \, 0 ) = X_2(p^2, \, p^2, \, 0) \ . \end{aligned}$$

(28)

The second condition is already built in the ansatz. Further, the ansatz (25)–(27) should be compatible with the scaling laws Eqs. (22) and (24). Recall that the gluon propagator $D(q^2)$ scales as $1/q^2$ at large $q^2$ and, therefore, $Y_1 (q^2)$ and $Y_3 (q^2)$ should approach a constant or scale, at large $q^2$, with a negative power of $q^2$ for large gluon momentum.

The explicit expressions for the longitudinal quark-gluon vertex form factors are derived from Eqs. (14), (25), (26) and (27). We can write the longitudinal components of the form factors as

$$\begin{aligned}&\lambda _1 = \frac{F(q^2)}{2} \Bigg \{ 2 \, B(p^2) \, D\left( \frac{p^2 + k^2}{2} \right) ~ Y_1(q^2) \nonumber \\&\quad +A(p^2) \left[ X_0(q^2) + \left( p^2 - p \cdot k \right) D\left( \frac{p^2 + k^2}{2} \right) ~ Y_3(q^2) \right] \Bigg \} \nonumber \\&\quad + (p \leftrightarrow k) \ , \nonumber \\&\lambda _2 = \frac{F(q^2)}{2 (k^2 - p^2)} A(p^2)\! \Bigg \{\left( p^2 + p\! \cdot k \right) D\left( \frac{p^2 + k^2}{2} \right) ~ Y_3(q^2) \nonumber \\&\quad - X_0(q^2) \Bigg \} + (p\leftrightarrow k) \ , \nonumber \\&\lambda _3 = \frac{F(q^2)}{p^2 - k^2} \Bigg \{ A(p^2) \left[ p^2 + p \cdot k \ \right] D\left( \frac{p^2 + k^2}{2} \right) ~ Y_1(q^2) \nonumber \\&\quad + B(p^2) \, X_0(q^2) \Bigg \}+ (p\leftrightarrow k), \nonumber \\&\lambda _4 = - \frac{F(q^2)}{2} D\left( \frac{p^2 + k^2}{2} \right) \Bigg \{ A(p^2) ~ Y_1(q^2) + B(p^2) ~ Y_3(q^2) \Bigg \} \nonumber \\&\quad - (p\leftrightarrow k) \, , \end{aligned}$$

(29)

where $A(p^2)=1/Z(p^2)$, $B(p^2)=M(p^2)/Z(p^2)$ with $Z(p^2)$ from Eq. (34) and $M(p^2)$ given in Eq. (35). The form factors $X_0 ( q^2)$, $Y_1(q^2)$ and $Y_3(q^2)$ are given later by Eq. (31). As written explicitly in Eqs. (29), all the symmetries coming from charge conjugation invariance [6] written in (10) are satisfied by our longitudinal form factors ansatz.

The use of the above ansatz with the gluon propagator as given in Eq. (32), of the ghost propagator as in Eq. (33) and of the quark functions given in (34) and (35) together with the longitudinal quark-gluon vertex (see Sect.  3.2) and the exact solution (14) of the Slavnov–Taylor identity (11) allows for the computation of a solution of the Schwinger–Dyson Equations, whose results will be discussed in the following sections.

4 Incorporating the Soft-Gluon Limit

To take into account the full QCD simulations results for the soft-gluon limit of $\lambda _1$ reported in [18] we insert (25)–(27) into (15) obtaining, after rotation to the Euclidean,

$$\begin{aligned} \lambda _1 (p^2) = \frac{F(0)}{Z(p^2)} \bigg \{ 1 ~ + ~ 2 \, M(p^2) \, D(p^2) \, Y_1(0) ~ - ~ 2 \, p^2 \, D(p^2) \, Y_3 (0) \bigg \} \ . \end{aligned}$$

(30)

A correlated fit of (30) to the lattice data using the ensemble with a pion mass of 295 MeV and a $\beta = 5.29$ [18], relying on the expressions for $D(p^2)$ and $M(p^2)$ given in [4], see also App. A, returns $Y_1(0) = 0.1726 \pm 0.0074$ GeV, $Y_3 (0) = -0.0806 \pm 0.008944$ and $F(0) / Z(0) = 1.211 \pm 0.029$ for a $\chi ^2/d.o.f. = 0.15$. In Fig.  1 we show both the lattice data and the fitted function (30). The above numbers for $Y_1(0)$ and $Y_3(0)$ provide a normalisation at zero momentum, that depends on the renormalisation scale $\mu$, for $Y_1(q^2)$ and $Y_3(q^2)$.

Fig. 1

Lattice data for $\lambda _1$ from full QCD simulations [18], in lattice units, and the fit to (30)

For the comparison of the fitted value for the ratio F(0) /Z(0) with the expression given in [4] and reproduced in Appendix A one has to consider the normalisation factor that comes from the renormalisation of $\lambda _1$. The fit was performed using bare lattice data and a gluon and ghost propagators renormalised in MOM-scheme at $\mu = 3$ GeV. This should be corrected before doing any comparison. One can estimate this global normalisation factor for $\lambda _1$ from the plot by demanding that $\lambda _1 (\mu ^2) = 1$. The quark renormalisation constant can be computed after setting $Z(\mu ^2) = 1$. These renormalisation constants allows the rescaling of the fitted ratio F(0) /Z(0) and only then estimate the value of F(0) that corresponds to a renormalised ghost propagator $F( \mu ^2 ) = 1$. It is this value, for $\mu = 3$ GeV, that should be compared with the ghost fitting function at zero momentum considered in [4]. It turns out that the fitted ghost propagator function used in [4] and that computed from the fit to the soft-gluon limit are compatible within one standard deviation. This gives us further confidence that the ansatz (25)–(27) for quark-ghost kernel form factors is able to capture the essential of the QCD dynamics.

5 Inverting the Schwinger–Dyson equations

The SDE are a set of two coupled equations that depend, for our ansatz, on three form factors. Furthermore, we would like to take into account the normalisation for $X_0$ given in Eq. (28), the corresponding UV limit and the normalisations for $Y_1(0)$ and $Y_3(0)$ from the fit to the soft-gluon limit of the lattice data, after appropriate rescaling to comply with the renormalisation scale. We recall the reader that in [4] we used $\mu = 4.3$ GeV as renormalisation scale and it will also be used here to solve the SDE. All these constraints can be taken in the calculation if all the functions are parametrised by Padé approximants

$$\begin{aligned} X_0 ( q^2)= & {} \frac{ 1 + a_{02} \, q^2 + a_{04} \, q^4}{1 + b_{02} \, q^2 + a_{04} q^4} \ , \nonumber \\ Y_1(q^2)= & {} \frac{Y_1(0) + a_{12} \, q^2 + a_{14} \, q^4 + a_{16} \, q^6 + a_{18} \, q^8}{ 1 + b_{12} \, q^2 + b_{14} \, q^4 + b_{16} \, q^6 + b_{18} \, q^8} \ , \nonumber \\ Y_3(q^2)= & {} \frac{Y_3(0) + a_{32} \, q^2 + a_{34} \, q^4 + a_{36} \, q^6 + a_{38} \, q^8}{ 1 + b_{32} \, q^2 + b_{34} \, q^4 + b_{36} \, q^6 + b_{38} \, q^8} \ . \nonumber \\ \end{aligned}$$

(31)

The coefficients in (31) were computed relying on simulating annealing to minimize the sum of the relative error of the scalar and vector equations. The numerical experiments show that it is relative easy to produce “solutions” whose maximum relative error for the SDE is of the order of 15%. However, for errors below the 10% value we found a single solution. As seen in Fig. 2 we found a solution that solves the SDE equations with a relative error, on each equation, below the 4% level. In the minimisation and to avoid poles on the Euclidean momenta real axis it was assumed that all the coefficients in the denominator are positive real numbers. We also report in Fig. 3 the quark wave function and the running quark mass, computed using our vertex solution, for the r.h.s. of the Schwinger–Dyson equations.

Fig. 2

Relative error for the Schwinger–Dyson equations for the solutions I and II reported in [4] for $\alpha _s = 0.22$ and the new solution considered here, computed using Padé approximants and taking into account the soft-gluon limit. In both solutions the propagators were renormalised at $\mu = 4.3$ GeV using the MOM-scheme

Fig. 3

The quark wave function and the running quark mass as given by the r.h.s. of the Schwinger–Dyson equations for our vertex solution, compared with the fit functions that parametrise $Z(p^2)$ and $M(p^2)$ and the lattice data for the quark propagator. For comparison the quark wave function was renormalised in the MOM-scheme at $\mu = 4.3$ GeV. The lattice data from $N_f = 2$ full QCD simulation in the Landau gauge [18, 22] for $\beta = 5.29$, $\kappa = 0.13632$ and for a $32^3 \times 64$ lattice was rescaled using the parameterisation reported in Appendix A. See [4] for the discussion on the parameterisation of the lattice quark propagator

Our parametrisation for $X_0(q^2)$ is the simplest Padé approximant that is compatible with the normalisation conditions $X_0(0) = X_0( + \infty ) = 1$ and allows for small deviations from unity as found in previous investigations [4, 10, 11]. Furthermore, taking as guide these previous calculations we expected a maximum of $X_0(q^2)$ below 1 GeV. Given that for small $q^2$, the function $X_0(q^2)$ is expected to grow, then $b_{02} < a_{02}$. If $X_0(q^2)$ has a maximum above 1 for $q < 1$ GeV, this demands $a_{02} < 1$ $\hbox {GeV}^{-2}$. All these constraints for $X_0$ were taking into account in the minimisation process.

In the minimisation of the error we also changed the powers of the numerator and denominator in the Padé approximants for $Y_1(q^2)$ and $Y_3(q^2)$ but only with those reported above we were able to find a solution of the SDE with a relative error below 4%. During the minimization process we observed that the first function to stabilize was $Y_1(q^2)$, followed by $Y_0(q^2)$ and then by $Y_3(q^2)$.

In Fig. 2 we show the relative error for the solution of the Schwinger–Dyson equations based on Padé approximants and the solutions reported in [4] computed with $\alpha _s = 0.22$. In all cases the relative error is below 4%.

Figure 3 shows $Z(p^2)$ and $M(p^2)$ computed from the r.h.s of the SDE when our vertex solution is used. The quark wave function follows very closely the original parameterisation used as input for the calculation, see App. A, with small deviations for momenta in the range $1 - 3$ GeV but reproducing the lattice data in the infrared region and the correct lattice data for larger momenta – see the discussion on the usage of the lattice quark propagator data in [4]. On the other hand, the running quark mass computed from the r.h.s. of the SDE is on top of the parameterisation for the same function used as input in the calculation.

The results reported in Fig. 3 help us to understand how the relative errors shown in Fig. 2 translate into the functions that describe the quark propagator. We recall that in Fig. 3 we are not solving the SDE for $Z(p^2)$ and $M(p^2)$ but are performing a consistency check of our approach and also on the ansatz (25)–(27) that we plug in into the solution (14) of the Slavnov–Taylor identity. Full expression for the quark-gluon form factors ansatz are given in (29).

6 Results

In Table 1 the coefficients for the solution that minimise the relative error of the SDE are reported. The corresponding form factors $X_0(q^2)$, $Y_1(q^2)$ and $Y_3(q^2)$ are shown in Fig. 4 and compared to the solutions computed in [4] with a completely different method, where the original SDE were replaced by Tikhonov regularised equations. All the represented solutions were computed using the same set of parameters, namely an UV hard cutoff of $\varLambda = 20$ GeV, $\alpha _s = 0.22$ and all propagators renormalised at $\mu = 4.3$ GeV. Moreover, for the various integrations, angular and momentum, we used exactly the same number of Gauss–Legendre points as in [4].

Table 1 Coefficients of the Padé approximant in (31) in powers of GeV. For each function, the upper line refers to the numerator coefficients in increasing power of $q^2$, while in the lower line are the coefficients for the denominator polynomial in increasing powers of $q^2$

Fig. 4

The solutions I and II found for $X_0(q^2)$, $Y_1(q^2)$ and $Y_3(q^2)$ in [4] and the new solution found with the Padé parameterisation. See text for details. Also shown are the 1-loop dressed perturbation theory for $X_0(q^2)$ (dashed lines) using both the tree level gluon-ghost vertex ($H(q^2)=1$) and an improved vertex ghost-gluon vertex [23]

For $X_0(q^2)$ the new solution is enhanced compared to those computed in [4], it has a maximum of $\sim$ 1.35 to be compared with $\sim$ 1.10 for the old solutions. The maximum of the new solution occurs at slightly larger $q \sim$ 450 MeV for the Padé based solution and $\sim$ 350 MeV for Tikhonov regularised solutions. The outcome of the one-loop dressed perturbation theory reported also in Fig. 4 have maxima that are similar to those of the Tikhonov regularised solution but occurring at a much larger scale, i.e. for $q \sim$ 1 GeV. The Padé based solution does not show any minima with $X_0(q^2) < 1$, as seen on the Tikhonov solutions, and approaches the UV normalisation condition $X_0(+ \infty ) = 1$ in a smoother way than the Tikhonov ones. In this respect the new solution follows closer the behaviour observed for the predictions of one-loop dressed perturbation theory.

The $Y_1(q^2)$ seen in Fig. 4 are quite similar up to $\sim 1$ GeV. The maximum of the Padé solution being slightly smaller than those of [4] and its deep infrared values, i.e. for $q \lesssim 200$ MeV, being larger to accommodate the lattice soft-gluon limit. For $q \gtrsim 1$ GeV, the various curves have similar structures, i.e. the same number of maxima and minima, but differ in UV. Here the Padé based solution approaches a negative constant value, while the Tikhonov solutions approach a positive constant value.

The Padé based solution for $Y_3(q^2)$ is different from those computed in [4]. It has a simplified structure that interpolates between its zero momentum value dictated by the lattice soft-gluon limit and a UV constant value that is about the same found for Sol. I in [4].

In [11] the authors solved simultaneously the SDE for the quark propagator together with the quark-ghost kernel, in its one-loop dressed perturbation theory formulation, to compute the various form factors $X_i$. At the qualitative level, but not quantitatively, our results point in the same direction. See [4] also for notation issues.

In Fig. 5 we report the various $\lambda _i$’s, the relevant quark-gluon vertex form factors, for the new Padé solution and compare them to the corresponding ones for Sol. I computed in [4]. Somehow surprisingly the differences between the two sets of form factors are minimal, with the exception of $\lambda _2$. These situations also occurs for other values of $\theta \ne 0$ (not represented here), the angle between the quark and the gluon momentum. This is a welcome feature, as the two solutions were computed in completely different and independent ways, giving confidence in our findings. In general the $\lambda _i$ based on the Padé solution for $X_0$ and $Y_i$ have slightly less structure as the oscillations observed in the Tikhonov regularised calculation are not present. $\lambda _1$ for the two solutions is very similar with the Padé based calculation showing a clear enhancement in the infrared region. The scale and the region of the absolute maximum of this form factor is the same for the two solutions. A similar conclusion holds also for $\lambda _3$. As Fig. 5 shows (see also the discussion in [4]) the main non-perturbative tensor structures of the quark-gluon vertex are associated with $\lambda _1$ and $\lambda _3$ and, in both cases, their maxima is at $p \sim \varLambda _{QCD}$. In what concerns the remaining form factors, the figure also shows that $\lambda _4$ of the Padé approximation is similar to that published for Sol. I in [4]. A large difference occurs for $\lambda _2$ that is infrared enhanced by about a factor of 2 and shows a single maximum instead of two for the [4] solution. For zero momentum, the value of $\lambda _1 (0)$ in the present calculation is about 2.1 times the one obtained for Sol. I of [4].

Our conclusion being that the inclusion of the information from the lattice soft-gluon limit in the calculation leads to an infrared enhancement of the quark-gluon vertex. The relative importante of the form factors found in [4] is also observed within the new solution. Our new calculation supports that the quark-gluon vertex is dominated by the form factors associated to the tree level vertex $\gamma _\mu$ and to $2 p_\mu + q_\mu$, with smaller contributions from the higher order tensor structures $(2 /p + /q) (2 p + q)_\mu$ and $\sigma _{\mu \alpha } (2 p + q)_\alpha$.

Fig. 5

Form factors for Sol. I reported in [4] (left) together with those associated with the new solution based on Padé approximantions (right)

To close our study, in Figs. 6, 7, 8 and 9 we illustrate the relative importance of the various quark-ghost kernel form factors to the quark-gluon $\lambda _i$ reported in Fig. 5 (right). For the dominant form factors $\lambda _1$ and $\lambda _3$ the major contribution for the infrared enhancement comes from $X_1$ ($\sim$ 72% and $\sim$ 83%), followed by $X_0$ ($\sim$ 20% and $\sim$ 15%) and a residual contribution from $X_3$ ($\sim$ 8% and $\sim$ 2%). The dominance of $X_1$ in the infrared region for these two form factors was unexpected. $\lambda _2$ is dominated by the $X_0$ contribution ($\sim$ 58%), followed by $X_3$ ($\sim$ 30%) and has a smaller contribution from $X_1$($\sim$ 12%). Finally, $\lambda _4$ is almost entirely built by the contribution coming from $X_1$. It is somehow unexpected that the various contributions to $\lambda _i$ do not follow the relative importance as suggested by perturbation theory. The terms proportional $X_1$ give the major contribution to $\lambda _3$, $\lambda _1$, $\lambda _4$ and $X_0$ is responsible for the major contribution only to $\lambda _2$.

Fig. 6

The contribution proportional to $X_0$ (left) and to $Y_1$ to the form factor $\lambda _1$

Fig. 7

The contribution proportional to $X_0$ (left) and to $Y_3$ to the form factor $\lambda _2$

Fig. 8

The contribution proportional to $X_0$ (left) and to $Y_1$ to the form factor $\lambda _3$

Fig. 9

The contribution proportional to $Y_1$ (left) and to $Y_3$ to the form factor $\lambda _1$

7 Discussion and summary

The lattice gluon propagator $D(p^2)$ appears as a multiplicative factor in the integral kernels of the Schwinger–Dyson Eqs. (16) and (17). If $D(p^2)$ changes, the vertex form factors need to change to compensate the modification of the propagator. For example, for a term as $D(p^2) \, \lambda _1$, if the propagator $D(p^2)$ is replaced by $D(p^2) f(p^2$), then $\lambda _1$ should change accordingly and $\lambda _1 \rightarrow \lambda _1/f(p^2)$. As there are indications that in full QCD $D(p^2)$ is IR suppressed [25, 26], meaning that $f(p^2)$ should be smaller than 1 in the infrared region, then the full QCD form factor $\lambda _1$ should increase by a factor $1/f(p^2)$ to keep the solution of the SDE for the self-energies unchanged. On the other hand, the gluon propagator also appears in the ansatzë for the quark-ghost kernel functions (25)–(27). This functional form is tied to the soft gluon limit of the quark-gluon vertex and the products appearing in (26) and (27) are fixed by the soft gluon limit. Therefore, one expects no or little modifications on the ansatz for the quark-ghost kernel. This seems to be consistent with the perturbative analysis of the quark-ghost kernel as the inclusion of the quark contributions are associated with higher order diagrams and are expected to be small. Hopefully this reasoning can be extended to lower momenta. This analysis of the contribution of the gluon propagator to the general solution of the SDE suggests that the form factors given in Fig. 5 underestimate the $\lambda _i$’s.

Our investigation of the solution of the Schwinger–Dyson equation takes into consideration only the longitudinal part of the quark-gluon vertex. Of course, the transverse part of the vertex also contributes to the dynamics of QCD and, in principle, its contribution to the quark gap equation is expected to be within the level of the errors found in the inversion performed herein. The transverse form factors have been discussed in the literature, see [27,28,29,30,31,32,33,34,35] and references therein, and, typically, they are written in a way to take into account the constraints coming from multiplicative renormalizability and gauge invariance.

In our description of the longitudinal form factors, see Eq. (29), setting $X_0=F(q^2)=1$, $Y_2=Y_3=0$ the vertex reproduces the longitudinal component of the Curtis-Pennington (CP) fermion-photon vertex [28]. The CP vertex satisfies the Ward-Green-Takahashi identity and is compatible with multiplicative renormalizability of the fermion SDE. The CP vertex also includes a transverse part, associated with the form factor $\tau _6$. The CP model was extended later taking into account the photon propagator SDE leading to the Kizilersü–Pennington (KP) vertex [30], that requires in addition to the CP vertex the inclusion of the transverse form factors $\tau _2$, $\tau _3$ and $\tau _8$. Our calculation shows that $\lambda _1$, $\lambda _3$ and $\lambda _4$ are dominated by $X_1$ against the expectation from the CP and KP vertex models, that suggests a dominance of $X_0$ contribution to the quark-ghost kernel. However, our results show, indeed, that the main contribution to $\lambda _2$ comes from $X_0$, in good agreement with the expectations coming from the CP and KP vertex models.

The anomalous chromomagnetic form factor, a contribution associated with the transverse form factor $\tau _5$ defined in Eqs. (8) and (9), was estimated in [31]. In the infrared region, $\tau _5$ attains a maximum of about 0.20 GeV$^{-1}$ and drops one order the magnitude for momentum of $\sim 1$ GeV. This form factor can be compared with the longitudinal tensor form factor $\lambda _4$ that, in the infrared region (below 1 GeV), takes values between − 1 and 1.5 GeV$^{-1}$. Of the four longitudinal form factors computed in the current work, the quark-gluon vertex is dominated by the contributions of $\lambda _1$ and $\lambda _3$, with $\lambda _4$ being subleading. However, our estimation of $\lambda _4$ is about 5 times larger than the anomalous chromomagnetic form factor obtained in [31]. If one takes this calculation as an estimation of the typical contribution of the transverse vertex to the quark gap equation, one expects small corrections that are due to the transverse form factors that, certainly, will decrease the longitudinal form factors obtained by inversion of the SDE.

Another interesting comparison can be done with the results for the quark-gluon vertex obtained in [34], were it is presented a self-consistent solution of the SDE for the quark propagator and the quark-gluon vertex function equation in a truncated form, that include only the three-gluon vertex and the dressed gluon in the decoupling scenario [36]. Although in their approach only the dependence on the gluon momentum is kept, their $\lambda _1$ form factor is comparable in magnitude to ours. In the infrared region, it attains the maximum value of about 3.5 in the chiral limit for the coupling constant having the value of $\alpha (\mu )=0.7427$ [34].The product of these two quantities is $\sim 2.6$ and can be compared to $\alpha _s \, \lambda _1 \, \sim 5$ at the peak of Fig. 5 with the two numbers being in the same ballpark. Recall that in our case the vertex depends on the gluon momentum, as in [34], but also on the quark momentum and on the angle between the quark and gluon momentum and, therefore, the comparison should be taken with care. If for our results one averages over the quark momentum, the estimation $\alpha _s \, \lambda _1$ will drop and the two results become closer. The differences in $Z(p^2)$ at the small-$p^2$ between our quark wave function, see Fig.  3, and that computed in [34] are overcome by the large vertex enhancements. Moreover, the lattice gluon dressing function is magnified in the IR with respect to the decoupling scenario used in the later work.

Our description of the quark wave function and of the the running mass reported in Fig. 3 should also be compared to the earlier studies of the quark propagator within the rainbow-ladder approximation. In [37, 38] a coupled set of renormalized Landau gauge truncated DSE’s for the quark, gluon, and ghost propagators was solved using a bare quark-gluon vertex, the Ball-Chiu and Curtis-Pennington models with and without nonabelian correction for the vertex. The framework used has a number of parameters that can be varied and it is possible to explore the parametric freedom of the models to obtain a running quark mass and a quark wave function that is similar to what is seen in Fig. 3. This is linked with a strong infrared enhancement of the kernel of the SDE. Our approach uses a clear grow of the gluon and ghost dressing functions, taken from recent lattice simulations, that via the Slavnov–Taylor identity and the modeling of the longitudinal form factors are fundamental to give the right strength of the SDE kernel in the infrared region.

The rainbow-ladder approach of Maris and Roberts [39], further developed by Maris and Tandy [40] proposes a quark–quark scattering kernel ansatzë, that has a strong IR enhancement to break dynamically the chiral symmetry, in a framework where the the axial-vector Ward-Takahashi identity and the Gell–Mann–Oakes–Renner relation are exactly satisfied. In addition the model preserves the one-loop renormalization group structure of QCD. This approach was explored in [41] and a fair description of the Lattice results for the running quark mass was achieved. Interesting enough, the study concluded that the dressed-gluon is not enough to obtain the existing lattice quark dressing data, unless the SDE kernel is enhanced for infrared momenta through the dressing of the quark-gluon vertex. Furthermore, the vertex in these works would correspond to the contribution of $\lambda _1$, which we found to have a quite strong IR enhancement as seen in Fig. 5; recall that in our case the form factors are function not only on the gluon momentum, but also on the quark momentum.

The symmetry-preserving Maris-Tandy model has been very successful in describing meson and baryon properties (see e.g. the reviews [42, 43]), emphasizing the role of the IR enhancement of the quark-gluon vertex as necessary to reproduce the hadron phenomenology.

In summary, the inclusion of the soft-gluon limit information in our analysis enhances the infrared region of the longitudinal form factors of the quark-gluon vertex. Furthermore, the relative contribution of the quark-ghost kernel form factors to the quark-gluon vertex do not follow in the infrared region their relative importance found in the ultraviolet region. Our solution for the longitudinal vertex form factors given in Eq. (29) together with the parametrization (31) and the results reported in Table  1, describes analytically all the components of the vertex considered, allowing easily further applications. It still remains a challenge to include the transverse components of the quark-gluon vertex. The error observed in the SDE, at the level of 4%, is expected to be associated not only with the limitation of the ansatz but also with the missing transverse form factors. The phenomenological implications of the present findings to the meson structure is unexplored and will be considered in a future publication.

Lattice ghost, quark and gluon propagators

For completeness, here we report on the fitting expressions used to model the lattice data for the various propagators. Note that, in all cases the model function incorporates the perturbative behaviour at high momenta. Details about how the model functions were computed and how well they reproduce the lattice propagator data can be found in [4].

The lattice Landau gauge gluon propagator and global fits to the data was published in [24]. In our analysis we use the following model function for the gluon propagador (see [24] for further details)

$$\begin{aligned} D (p^2) = Z \, \frac{ p^2 + M^2_1}{p^4+ M^2_2 \, p^2 + M^4_3} \left[ \omega \, \ln \left( \frac{p^2 + m^2_0}{\varLambda ^2_{QCD}}\right) + 1 \right] ^{\, \gamma } \ , \nonumber \\ \end{aligned}$$

(32)

where the gluon anomalous dimension is $\gamma = -13/22$, $Z = 1.36486 \pm 0.00097$, $M^2_1 = 2.510 \pm 0.030 \hbox { GeV}^2$, $M^2_2 = 0.471 \pm 0.014 \hbox { GeV}^2$, $M^4_3 = 0.3621 \pm 0.0038 \hbox { GeV}^4$, $m^2_0 = 0.216 \pm 0.026 \hbox { GeV}^2$. In Eq. (32) $\omega = 33 \, \alpha _s ( \mu ) / 12 \pi$ and to reproduce the lattice data $\varLambda _{QCD} = 0.425 \hbox { GeV}$ and the strong coupling constant reads $\alpha _s ( \mu = 3 \hbox { GeV} ) = 0.3837$. We call the reader attention that in the inversion of Dyson-Schwinger equation the propagator is renormalised, in the MOM scheme, at a difference scale. This implies the computation of a global renormalisation constant that multiplies the above $D(p^2)$. The same reasoning applies to the ghost propagator and to the quark wave function.

For the ghost propagator we take the data reported in [21] for the simulation using a $80^4$ and use the following functional form to describe the lattice data

$$\begin{aligned}&D_{gh} (p^2) = \frac{F(p^2)}{p^2} \nonumber \\&\quad = \frac{Z}{p^2} ~ \frac{ p^4 + M^2_2 \, p^2 + M^4_1}{p^4+ M^2_4 \, p^2 + M^4_3} \times \nonumber \\&\qquad \times \, \left[ \omega \, \ln \left( \frac{p^2 + \frac{m^4_1}{p^2 + m^2_0}}{\varLambda ^2_{QCD}}\right) + 1 \right] ^{\, \gamma _{gh}} . \end{aligned}$$

(33)

The fit to the lattice data results in $Z = 1.0429 \pm 0.0054$, $M^4_1 = 18.2 \pm 5.7 \hbox { GeV}^4$, $M^2_2 = 33.4 \pm 6.4 \hbox { GeV}^2$, $M^4_3 = 6.0 \pm 2.7 \hbox { GeV}^4$, $M^2_4 = 29.5 \pm 5.7$ $\hbox {GeV}^2$, $m^4_1 = 0.237 \pm 0.049$, $m^2_0 = 0.09 \pm 0.42 \hbox { GeV}^2$. The ghost anomalous dimension reads $\gamma _{gh} = - 9/44$ and $\omega$ and $\varLambda _{QCD}$ take the same values as in the gluon fitting function (32).

For the quark propagator we consider the result of a $N_f = 2$ full QCD simulation in the Landau gauge [18, 22] for $\beta = 5.29$, $\kappa = 0.13632$ and for a $32^3 \times 64$ lattice. For this particular lattice setup, the corresponding bare quark mass is 8 MeV and the pion mass reads $M_\pi = 295$ MeV. For the quark wave function $Z(p^2)$ we corrected the lattice data to become compatible with perturbation theory at high momenta, see [4] for details, and use the following functional form

$$\begin{aligned} Z(p^2) = Z_0 ~ \frac{ p^4 + M^2_2 \, p^2 + M^4_1 }{p^4 + M^2_4 \, p^2 + M^4_3} \end{aligned}$$

(34)

with $Z_0 = 1.11824 \pm 0.00036$, $M^4_1 = 1.41 \pm 0.18 \hbox { GeV}^4$, $M^2_2 = 6.28 \pm 1.00 \hbox { GeV}^2$, $M^4_3 = 2.11 \pm 0.28 \hbox { GeV}^4$, $M^2_4 = 6.20$ $\pm 0.98$ $\hbox {GeV}^2$. The choice of the functional form for the running quark mass was

$$\begin{aligned} M(p^2) = \frac{ m_q ( p^2) }{ \left[ A + \log ( p^2 + \lambda \, m^2_q(p^2) ) \right] ^{\gamma _m} } \end{aligned}$$

(35)

where $\gamma _m = 12/29$ is the quark anomalous dimension for $N_f = 2$ and

$$\begin{aligned} m_q(p^2) = M_q ~\frac{ p^2 + m^2_1}{ p^4 + m^2_2 \, p^2 + m^4_3 } + m_0 \ . \end{aligned}$$

(36)

The parameters in Eqs. (35) and (36) are $M_q = 349 \pm 10 \hbox { MeV GeV}^2$, $m^2_1 = 1.09 \pm 0.43 \hbox { GeV}^2$, $m^2_2 = 0.92 \pm 0.28 \hbox { GeV}^2$, $m^4_3 = 0.42 \pm 0.15 \hbox { GeV}^4$, $m_0 = 10.34 \pm 0.63$ MeV and $A = -2.98 \pm 0.25$ and $\lambda = 1 \hbox { GeV}^2/\hbox {MeV}^2$.