The Quark-Gluon Vertex and the QCD Infrared Dynamics

The Dyson-Schwinger quark equation is solved for the quark-gluon vertex using the most recent lattice data available in the Landau gauge for the quark, gluon and ghost propagators, the full set of longitudinal tensor structures in the Ball-Chiu vertex, taking into account a recently derived normalisation for a quark-ghost kernel form factors and the gluon contribution for the tree level quark-gluon vertex identified on a recent study of the lattice soft gluon limit. A solution for the inverse problem is computed after the Tikhonov linear regularisation of the integral equation, that implies solving a modified Dyson-Schwinger equation. We get longitudinal form factors that are strongly enhanced at the infrared region, deviate significantly from the tree level results for quark and gluon momentum below 2 GeV and at higher momentum approach their perturbative values. The computed quark-gluon vertex favours kinematical configurations where the quark momentum $p$ and the gluon momentum $q$ are small and parallel. Further, the quark-gluon vertex is dominated by the form factors associated to the tree level vertex $\gamma_\mu$ and to the operator $2 \, p_\mu + q_\mu$. The higher rank tensor structures provide small contributions to the vertex.


Introduction
The interaction of quarks and gluons is described by Quantum Chromodynamics [1][2][3][4], a renormalizable gauge theory associated to the color gauge group SU (3). The quarkgluon vertex has a fundamental role in the description of hadron phenomenology, to the understanding of chiral symmetry breaking mechanism and, eventually, on the realisation of confinement. Despite its relevance for strong interactions, our knowledge of the quarkgluon vertex from first principles calculations is relatively poor. At the perturbative level, only recently a full calculation of the twelve form factors associated to this vertex was published [5]. The twelve form factors needed to describe the quark-gluon vertex have been computed for various kinematical configurations namely the symmetric configuration (equal incoming, outgoing quark and gluon squared momenta), the on-shell configuration (quarks on-shell with vanishing gluon momentum) and the asymptotic limit (all momenta much larger than the current quark mass). The study of the asymptotic limit of the vertex was then used to investigate various ansätze that can be found in the literature [6][7][8][9][10][11] and, in particular, test its description of the ultraviolet regime. At the non-perturbative level, the quark-gluon vertex has been studied within continuum approaches to QCD by several authors [11][12][13][14][15][16][17][18][19][20][21][22][23]. Typically, the computation is performed after writing the vertex in terms of other QCD vertices and propagators, preferably taking into account the perturbative tail, in order to simplify its calculation. Most of the computations also include only a fraction of the total number of the twelve form factors required to describe the quark-gluon vertex. In the calculation performed within massive QCD, i.e. using the Curci-Ferrari model [18], the vertex was computed in perturbation theory and all its (perturbative) tensor structures were accessed. In [17,22,23] the authors solved the theory at the non-perturbative level gathering information on the vertex from QCD symmetries and relying on one-loop dressed perturbation theory.
The quark-gluon vertex has also been computed using lattice simulations both for quenched QCD [24][25][26] and for full QCD [27]. On lattice simulations typically only a limited set of kinematical configurations are accessed with the soft gluon limit, defined by a vanishing gluon momenta, being the most explored. In particular, for full QCD so far only a single form factor, that associated with the tree level tensor structure, was measured in the soft gluon limit. The lattice calculations need a proper estimation of the lattice artefacts, so that they can be subtracted to the results of the simulations, to report the corresponding continuum functions.
One can also find on the literature attemps to combine continuum non-perturbative QCD equations with results from lattice simulations to study the quark-gluon vertex. In [28] a generalized Ball-Chiu vertex was used in the quark gap equation together with lattice results for the quark, gluon and ghost propagators to address the quark-gluon vertex. In [29], the full QCD lattice data for λ 1 was investigated relying on continuum information about the vertex. The use of continuum equations with results coming from lattice simulations requires high quality lattice data to feed the continuum equations that should be solved for the vertex. In this approach, the computation of a solution of the continuum equation requires assuming some type of functional dependence for various functions. In recent years, there has been an effort to improve the quality of the lattice data, in the sense of being closer to the continuum and producing simulations with large statistical ensembles, both for propagators and for vertex functions.
For the practitioner oftentimes it is sufficient to have a good model of the vertex that should incorporate the perturbative tail to describe the ultraviolet regime, some "guessing" for the infrared region and, hopefully, comply with perturbative renormalization [5][6][7][8]31]. A popular and quite successful model was set in [32], named the Maris-Tandy model. It assumes a unique tensor basis for the vertex, its tree level tensor structure, and simplifies its functional form. Further, this model assumes that the vertex depends only on the gluon momentum and accommodates the results from perturbative theory at high momentum transfer with a considerable enhancement of the vertex at low momenta. Such type of vertex that appears in the Dyson-Schwinger and the Bethe-Salpeter equations can be seen as a reinterpretation of the full vertex tensor structure, after rewriting its main components in a way that formally can be associated with the gluon propagator. Although the Maris-Tandy model is quite successful for phenomenology, the model is not able to describe the full hadronic properties and fails to explain the mass splittings of the ρ and a 1 parity partners, underestimates the weak decay constants of heavy-light mesons and cannot reproduce simultaneously the mass spectrum and decay constants of radially excited vector mesons to point out some known limitations. For a more complete description see, for example, [33,34] and references therein.
The goal of the present work is to explore further the quark-gluon vertex in the nonperturbative regime from first principles calculations. Our approach follows the spirit of the calculations initiated in [28] that combines continuum methods with results from lattice simulations to solve the quark Dyson-Schwinger equation for the vertex. If in that work the quark-gluon vertex was computed using a generalised Ball-Chiu vertex, herein we go much beyond and include other tensor structures that imply taking into account the dependence of the vertex with the quark momentum and the angle between quark and gluon momentum. As in the above cited work, the calculations performed use the Landau gauge to profit from the recent high quality lattice data for the quark, the gluon and the ghost propagators and taking into consideration only the longitudinal components of the vertex. Moreover, the computation also incorporate the recent analysis of the full QCD lattice simulation for the quark-gluon vertex in the soft gluon limit that identifies an important contribution associated to the gluon propagator [29]. We use a Slavnov-Taylor identity to write the longitudinal components of the vertex as a function of the quark wave function, the running quark mass, the quark-ghost kernel form factors and the ghost propagator. We also incorporate the normalisation of the quark-ghost kernel form factors X 0 , see below for definitions, derived in [17] for the soft gluon limit. The normalisation of X 0 had an important role in the analysis of the full QCD lattice data analysis for λ 1 , the form factor associated with the tree level tensor structure (see below for definitions), done in [29].
Our solution for the quark-gluon vertex returns a X 0 that deviates only slightly from the normalisation condition referred above and longitudinal quark-gluon vertex form factors that are strongly enhanced in the infrared region. The enhancement of the four longitudinal form factors occurs for quark and gluon momentum below 2 GeV and at high momentum the form factors approach their perturbative values. The computed quark-gluon vertex is also a function of the angle between the quark four momentum p and the gluon four momentum q that, clearly, favours kinematical configurations where p and q of the order of 1 GeV or below. Furthermore, we found that the vertex is enhanced when all momenta entering the vertex (see Fig. 1) tends to be parallel in pairs, which is solved by the compromisse that the momenta are restricted to a region around Λ QCD . Within our solution for the quark-gluon vertex the dominant form factors are those associated to the tree level vertex γ µ and to the scalar 2 p µ + q µ , with the higher rank tensor structures giving subleading contributions.
The paper is organised as follows. In Sec. 2 we introduce the notation for the propagators, the Dyson-Schwinger equations and the quark-gluon vertex. Moreover, we use a Slavnov-Taylor identity to rewrite the vertex in terms of the quark propagators functions Figure 1: The quark-gluon vertex. and the quark-ghost kernel. The parametrization of the quark-ghost kernel is also discussed. In Sec. 3 the scalar and vector components of the DSE in Minkowsky space are given together with the corresponding kernels and we start to set up the ansatz to be used to solve the integral equations. In Sec. 4 the DSE are rewritten in Euclidean space and, after performing a scaling analysis of the integral equations, we introduce the ansatz for the vertex. In Sec. 5 we give the details of the lattice data used in the current work for the various propagators and the details of the functions that parametrize the lattice data. The kernels for the Euclidean space DSE are discussed in Sec. 6, together with the solutions for the vertex of the gap equation. The quark-gluon vertex form factors are reported in Sec. 7 for several kinematical configurations. Finally, on Sec. 8 we summarise and conclude.

The Quark Gap Equation and the Quark-Gluon Vertex
In this section the notation used through out the article is defined. The equations to be discussed below and in the first part of this work have all expressions defined in Minkowsky space with the diagonal metric g = (1, −1, −1, −1). Following the notation of [35], in the quark-gluon vertex represented in Fig. 1 all momenta are incoming and, therefore, verify (2.1) The one-particle irreducible Green's function associated to the vertex reads where g is the strong coupling constant and t a are the color matrices in the fundamental representation.
The quark propagator is diagonal in color and its spin-Lorentz structure is given by where Z(p 2 ) = 1/A(p 2 ) stands for the quark wave function and M (p 2 ) = B(p 2 )/A(p 2 ) is the renormalization group invariant running quark mass. The inverse quark propagator reads

4)
[ p Figure 2: The Dyson-Schwinger equation for the quark. The solid blobs denote dressed propagators and vertices.
The Dyson-Schwinger equation for the quark propagator, also named the quark gap equation, is represented in Fig. 2 and can be written as where Z 2 is the quark renormalization constant, m bm the bare current quark mass and the quark self-energy is given by where Z 1 is a combination of several renormalization constants, ∆ ab µν (q) is the gluon propagator which, in the Landau gauge, is given by In the following, both ∆ ab µν (q) or the form factor ∆(q 2 ) will be referred as the gluon propagator.
A key ingredient in gap equation (2.5) is the quark-gluon vertex. Indeed, it is only after knowing Γ a µ that one can compute Z(p 2 ) and M (p 2 ). The Lorentz structure of the quarkgluon vertex Γ µ , see Eq. (2.2), can be decomposed into longitudinal Γ (L) and transverse Γ (T ) components relative to the gluon momenta, i.e. one writes where, by definition, By choosing a suitable tensor basis in the spinor-Lorentz space, Γ µ can be written as a sum of scalar form factors that multiply each of the elements of the basis. The full vertex Γ µ requires twelve independent form factors and, using the Ball and Chiu basis [6] it becomes The operators associated to the longitudinal vertex are while those associated to the transverse part of the vertex read

QCD Symmetries and the Quark-Gluon Vertex
The global and local symmetries of QCD constrain the full vertex Γ µ and connect several of the Green's functions theory. For example, the global symmetries of QCD require that the form factors λ i and τ i to be either symmetric or anti-symmetric under exchange of the two first momenta -see, e.g., ref. [35] and references therein. On the other hand, gauge symmetry imply that the Green functions also satisfy the Slavnov-Taylor identities (STI) [36,37]. This identities play a major role in our understanding of QCD and, in particular, the longitudinal part of the quark-gluon vertex is constrained by the following identity where the ghost-dressing function F (q 2 ) is related to the ghost two-point correlation function as and H and H are associated to the so-called quark-ghost kernel. As discussed in [35], these functions can be parametrized in terms of four form factors as The STI given in Eq. (2.14) can be solved with respect to the vertex [13] and, in this way, the longitudinal form factors λ i (p 1 , p 2 , p 3 ) can be written in terms of the quark propagator functions A(p 2 ), B(p 2 ) and the quark-ghost kernel functions X i and X i as A nice feature of the above solution for the various form factors that can be check by direct inspection is that the symmetry requirements on the λ i due to charge conjugation are automatically satisfied independently of the functions A, B, X i and X i . This is a particularly important point if one aims to model the vertex.

Decomposing the Dyson-Schwinger Equation into its Scalar and Vector Components
The Dyson-Schwinger equation for the quark propagator is written in (2.5), with the quark self-energy being given by (2.6). This equation can be projected into its scalar and vector components by taking appropriate traces. The scalar part of the equation is given by the trace of (2.5) which, after some algebra, reduces to after insertion of the vertex decomposition given in (2.8), taking into account only its longitudinal part, where k = p − q, and C F = 4/3 is the Casimir invariant associated to the SU(3) fundamental representation. The vector component of (2.5) is obtained after multiplication by / p and then taking the trace of the resulting equation to arrive on The two equations (3.1) and (3.3) can be simplified further by modelling the quarkgluon vertex. For example, in [13,28] the vertex was parametrized using the solution of the Slavnov-Taylor identity (2.17)-(2.20) and setting X 1 = X 2 = X 3 = 0. The rationale for such a choice comes from perturbation theory which gives, at tree level, X 0 = 1 and X 1 = X 2 = X 3 = 0. This ansatz, that ignores all form factor associated to the quarkghost kernel but X 0 , assumes that at the non-perturbative level the hierarchy of the form factors follows its relative importance that is observed for high momentum. Furthermore, in order to compute a solution of the Dyson-Schwinger equations it was assumed a further restriction on X 0 , namely that it depends only on the incoming gluon momenta, i.e. that X 0 = X 0 (q 2 ).
In order to proceed the analysis of the Dyson-Schwinger equations for the vertex it will be assumed that, in what concerns the momentum dependence, the form factors associated to the quark-ghost kernel factorize as where g i (p 2 1 , p 2 2 ) = g i (p 2 2 , p 2 1 ) are symmetric functions of its arguments. This type of factorisation is compatible, for example, with the Maris-Tandy quark-gluon vertex [32] and simplifies considerably the analysis of the solutions of the Dyson-Schwinger equations.
In [17] it was proved that, to all-orders, that, for the ansatz considered above, imply g(p 2 , p 2 )X 0 (0) = 1 and g 1 (p 2 , p 2 )X 1 (0) = g 2 (p 2 , p 2 )X 2 (0) . (3.6) In order to comply with the second relation given in Eq (3.5) it will be assumed from now on that X 1 = X 2 for any kinematical configuration. Note also that by choosing the g i to be symmetric functions of the arguments, the form factors X i and X i become identical. By taking into account the ansatz for the quark-gluon vertex just set in into the solutions of the Slavnov-Taylor identities (2.17)-(2.20), it follows that where and k = p − q; we call the reader's attention that we have used X i to representX i . Then, the scalar component of the Dyson-Schwinger equations becomes where the kernels are defined as with the kernels given by

The Dyson-Schwinger Equations in Euclidean Space
Our aim is to solve the Dyson-Schwinger equations for the quark-ghost kernel, said otherwise for the quark-gluon vertex, and this requires the knowledge of the quark, gluon and ghost propagators. For that we will rely on lattice inputs that provide first principles nonperturbative results and also demand that the above expressions should be rewritten in Euclidean space. The Wick rotation to go from Minkowski to Euclidean space is achieved by making use of the following substitutions on Eqs. (3.12) and (3.16). For completeness, below we will provide all expressions in Euclidean space.
The scalar component of the Dyson-Schwinger equations reads while its vector component is given by The kernels appearing in Eqs. (4.2) and (4.3) are The study of the solutions of the above equations using lattice inputs requires the use of the renormalized Dyson-Schwinger equations and, therefore, all quantities appearing on these equations should be finite. This requirement constrains the integrand functions g i (p 2 , (p − q) 2 ) X i (q 2 ) and, in particular, its possible behaviour in the limits where q → 0 and p → +∞.
Let us start by considering the ultraviolet limit of the integrand functions appearing in Eqs. (4.2) and (4.3). In the large q limit it follows that up to logarithmic corrections associated to the various propagators. In this limit, the integrand function appearing on the scalar equation (4.2) read The requirement of having a finite integral demands that at large q or that these functions are proportional to a higher negative power of q. The logarithmic corrections, not taken into account in this analysis, are sufficient to avoid the UV logarithmic divergence suggested by the naive power counting. Indeed, these logarithmic corrections that are introduced by the renormalization group improvement are, for large momenta, of type log(q 2 /Λ 2 ) γ , with γ being one of the anomalous dimensions. Our large q analysis should take into account the logarithmic corrections coming from the gluon, the ghost and the quark propagators that for N f = 2 result in γ = γ glue + γ ghost + γ quark = −137/116. Then, assuming a large q behaviour as in (4.12) times the log correction, the integration function at high momenta becomes resulting in a finite value for the integral. The difference between the naive power counting and taking into account the log corrections is illustrated on Fig. 3. As seen on this figure, the log corrections suppress further the integrand function at high momenta.
In what concerns the quark-ghost kernel form factor X 0 (q 2 ) at high energies, the power counting analysis is compatible with a X 0 (q 2 ) = 1 at large momenta as required by perturbation theory and by the all-orders result summarised in Eq. (3.6).
The same analysis for the vector component (4.3) gives, up to logarithmic corrections, where {· · · } stand for finite expressions involving A(p 2 ), A(q 2 ), B(p 2 ) and B(q 2 ). The conditions given in (4.12) are sufficient to ensure a finite result associate to the UV integration over q on the vector component of the Dyson-Schwinger equations. The QCD dynamics generates infrared mass scales for the quark and gluon propagators that are sufficient to eliminate possible infinities associated to the low momentum limit in the integral of the quark gap equation.
For full QCD, the λ 1 form factor of the quark-gluon vertex was computed in the soft gluon limit, i.e. vanishing gluon momenta, using lattice simulations in [27]. The analysis of the lattice data performed in [29] show that the lattice data is well described by where a and b are constants, that in terms of X 1 and X 3 translates into . This results suggests to write that in the high q limit gives X 1 ∼X 1 (q 2 )/q 2 and regularizes the ultraviolet behaviour in agreement with the discussion summarised in (4.12). Similarly, equation (4.16) also suggests to write and at high q momenta X 3 ∼X 3 (q 2 )/q 2 and, in this way, the ultraviolet problems referred in (4.12) are solved. Furthermore, for large quark momentum the ansatz (4.17) and (4.18) give implying the vanishing of the kernels (4.4) -(4.9) for sufficiently large p.
In summary, our ansatz for the quark-gluon vertex used to solve the Dyson-Schwinger equations reads (4.20) The quark gap equation should be solved taking into account the constraint (3.6) that demands X 0 (0) = X 0 (q → +∞) = 1 .  Our goal is to use the gap equation to study the quark-gluon vertex and, therefore, the knowledge of the various propagators over all range of momenta appearing in the integral equation is required. This is achieved by fitting the Landau gauge lattice propagators using model functions that are compatible with 1-loop renormalization group improved perturbation theory. In this way, it is ensured that the perturbative tails are taken into account properly in the parameterization of the various functions used to solve the Dyson-Schwinger equations. In order to compare the present work with the results of [28], on App. B we compared the new fits, to be discussed below with those previously used. Note the differences between the two sets of curves that necessarily change quantitatively, but not qualitatively, the solutions presented herein and in [28].

Landau gauge lattice gluon and ghost propagators
The lattice gluon propagator has been computed in the Landau gauge both for full QCD and for the pure Yang-Mills. Nowadays, the gluon propagator is well known for the pure Yang-Mills theory and it was calculated in [42] for large statistical ensembles and for large physical volumes ∼ (6.6 fm) 4 and ∼ (8.2 fm) 4 . Furthermore, in [42] the authors provide global fits to the lattice data that reproduce the 1-loop renormalisation group summation of the leading logarithmic behaviour. Of the various expressions, we will use to solve the Dyson-Schwinger equations the fit to the (6.6 fm) 4 volume given by the χ 2 /d.o.f. However, given that the level of precision achieved on lattice simulations for the quark propagator is considerable smaller than for the gluon propagator, one should distinguish between the various fitting functions provided in [42]. Our option was to use the simplest functional form given in this work. The lattice data for the Landau gauge gluon dressing function p 2 ∆(p 2 ), renormalized in the MOM-scheme at the mass scale µ = 3 GeV and the fit associated to Eq. (5.1) can be seen on the left part of Fig. 4.
For the ghost propagator we take the data reported in [41] for the 80 4 lattice simulation and fit the lattice data to the functional form

Lattice Quark Propagator
For the quark propagator we consider the N f = 2 full QCD simulation in the Landau gauge of [27] for β = 5.29, κ = 0.13632 and for a 32 3 × 64 lattice. For this particular lattice setup, the corresponding bare quark mass is 8 MeV and the pion mass reads M π = 295 MeV.
Our fittings to the lattice data, see below, take into account that the lattice data is not free of lattice artefacts; see [27] for further details. At high momenta the lattice quark wave function Z(p 2 ) is a decreasing function of momenta, a behaviour that is not compatible with perturbation theory that predicts a constant Z(p 2 ) in the Landau gauge.
As reported in [27], the analysis of the lattice artefacts relying on the H4 method suggests that, indeed, Z(p 2 ) is constant at high p. In order to be compatible with perturbation theory, we identify the region of momenta where Z(p 2 ) is constant and, for momenta above this plateaux, we replace the lattice estimates of Z(p 2 ) by constant values, i.e. the higher value of the quark wave function belonging to the plateaux. The original lattice data and the ultraviolet corrected lattice data can be seen on the left of Fig. 5. The UV corrected lattice data is then fitted to the rational function The removal of the lattice artefacts for the running quark mass is more delicate when compared to the evaluation of the quark wave function lattice artefacts [24,38]. The lattice data published in [27] and reported on Fig. 5 (right) was obtained using the so called hybrid corrections to reduce the lattice effects [24] . The hybrid method results in a smoother mass function when compared to the one obtained by applying the multiplicative corrections. The differences on the corrected running mass between the two methods occur for momenta above 1 GeV, with the multiplicative corrected running mass being larger than the corresponding hybrid estimation. The running mass provided by the two methods, corrected for the lattice artefacts, seems to converge to the same values at large momentum.
The running mass reported on Fig. 5 (right) is not smooth enough to be fitted. To model the lattice running mass in a way that reproduces the ultraviolet and the infrared lattice data and is compatible with the perturbative behaviour at high moment, we remove some of the lattice data at intermediate momenta. On Fig. 5 the data in the region with an orange background was not taken into account in the global fit to the running quark mass. The remaining lattice data was fitted to where γ m = 12/29 is the quark anomalous dimension for N f = 2 and The fitted parameters are M q = 349 ± 10 MeV GeV 2 , m 2  integration will be performed as described in App. A, i.e. by introduction an hard cutoff Λ, and with all integrations performed with Gauss-Legendre quadrature. For the angular integration we consider 500 Gauss-Legendre points as in [28]. After angular momentum integration, one is left with the kernels that can be seen on Figs. 6, 7 and 8 without taking into account the Gauss-Legendre weights associated to the integration over the gluon momentum. The inclusion of the Gauss-Legendre weights associated to the q momentum integration does not change the outcome reported on Figs. 6, 7 and 8 and the main difference being that the associated numerical values are considerable smaller. The major contributions of the N A reproduce the behaviour observed in [28] and, therefore, the momentum integration in Eqs. (4.2) and (4.3) associated to the kernels that are coupled to X 0 (q 2 ) is finite.
The function N A (p, q) displays a similar pattern and, again, the integration over the gluon momentum associated with N (1) A is expected to be well behaved. On the other hand the remaining kernels, i.e. N (1) , are all increasing functions of q. The requirement of a finite integration over q demands that X (1) and X (3) approach zero fast enough to compensate the increase with q of these kernel functions; see the discussion of the kernels ultraviolet limit in Sec. 4. The ansatz (4.20) -(4.22) adds a multiplicative gluon propagator term that is just enough to regularize the high momentum associated to N A (p, q), i.e. their main contribution to the integral equations is for p 2 GeV and q 2 GeV, and, therefore, the inclusion of the gluon propagator in the kernels makes the integration over q finite.
The Dyson-Schwinger equations are solved using a hard cutoff that is set to Λ = 20 GeV and all quantities are renormalized in the MOM scheme and using the same renormalization scale as in [28], i.e µ = 4.3 GeV, so that one can compare the results of the two works. It follows that the renormalized quantites satisfy the identities The bare quark mass quoted in the lattice simulation for the ensemble used here reads m bm = 8 MeV [27]. In the following we set Z 1 = 1, take the value for Z 2 from the vector component of the gap equation at the cutoff and "measure" the bare quark mass using the A (p, q)/p 2 kernels including the term of the gluon propagator as defined in (4.21). See also the caption of Fig. 6.
A (p, q)/p 2 kernels including the term of the gluon propagator as defined in (4.22). See also the caption of Fig. 6. scalar component of the gap equation at the cutoff momenta. In this way m b.m. does not coincide with the value quoted in the simulation but, as can be seen below, its value is close to the 8 MeV quoted above.
The results shown on Secs. 6.1, 6.2, 6.3 and 6.4 were computed using the same value for α s (µ) = 0.295 as in [28]. In Sec. 6.5 we allow α s (µ) to change and provide a best value. From Sec. 6.5 onwards, the results reported use the optimal value for the strong coupling constant.

One-Loop Dressed Perturbation Theory for
The ansatz considered here demands the knowledge of the three form factors X 0 , X 1 and X 3 . However, the quark gap equation provides only two independent equations and, therefore, it is not possible to compute all the form factors at once for the full range of momenta.
A first look at the quark-ghost kernel form factors is possible if one computes X 0 within one-loop dressed perturbation theory with a simplified version of a quark-ghost kernel where Figure 11: One-loop dressed perturbation theory one sets X 1 = X 3 = 0 and, then, solve the gap equation to estimate X 1 and X 3 . The way the solutions of the Dyson-Schwinger equations for X 1 and X 3 are built also illustrates the numerical procedure used to solve the integral equations.
The one-loop dressed approximation to the quark-ghost kernel is represented on Fig. 11 that, in the simplified version of kernel, translates into the following integral equation and where H 1 (q 2 ) stands for the ghost-gluon vertex. When solving this equation we consider two version of the ghost-gluon vertex, namely its tree level version where H 1 (q 2 ) = 1 and an enhanced dressed vertex as given in [43] that takes 3) to a linear system of equations that was solved using the QR decomposition of the matrix appearing in the linear system. The numerical solutions for X 0 can be seen on Fig. 12 and are, essentially, those reported in [28]. According to one-loop dressed perturbation theory, the deviations of X 0 (q 2 ) from its tree level value are, at most, of the order of 15%. We have also looked at the iterative solutions for (6.3). No convergence was observed and the solution produced by a single iteration results in a X 0 (q 2 ) that is enhanced relative to the solutions of Fig. 12.
The estimation of X 0 allows to solve the gap equation for X 1 and X 3 . In order to compute a solution of the Dyson-Schwinger equations, after performing the angular integration, the scalar and vector components of the equations are rewritten in the form of the larger  Figure 12: Simplified one-loop dressed perturbation theory estimation for X 0 (q 2 ) from now on we will adopt the short name version B = N X to refer to this linear system of equations. Note that X 1 has mass dimensions, while X 0 and X 3 are dimensionless. We could have rescaled X 1 to make it dimensionless. The only natural scale in QCD is provided by Λ QCD or, alternatively, by the propagators at some mass scale. However, given that we have no idea of which mass scale to consider, we did not introduced any scale or, said otherwise, the results reported here can be viewed as taking this mass scale to be 1 GeV. This is also true for the full solution of the Dyson-Schwinger equations discussed in next section.
A direct solution of B = N X results in a meaningless result, with the components of X oscillating over very large values due to the presence of very small eigenvalues of the matrix N , that translates the ill defined problem in our hands. The linear system can be solved using the Tikhonov regularization [44] that is equivalent to minimize ||B − N X|| 2 + ||X|| 2 , where is a small parameter to be determined in the inversion. In this way we look for solutions that solve the linear system but whose norm is small. For real symmetric matrices, Tikhonov regularization replaces the original system by N T B = (N T N + )X. In our case N is not a symmetric matrix and we will solve the system as given by its normal form.
In order to determine the regularization parameter, we solve N T B = (N T N + )X for various values of and look at how ||B − N X|| 2 and ||X|| 2 change with . The outcome of the inversions for different can be seen on Fig. 13. On this figure, smaller values of are closer to the original ill defined problem and correspond to solutions with larger norms for X 1 and X 3 . On the other hand, larger values of are associated to the solutions of the modified linear system with smaller X 1 and X 3 norms. The optimal value of is given by the solution whose residuum, i.e. the difference between the lhs and the rhs of the original  Figure 13: Residuum versus norm for the scalar and vector equation when solving the gap equation for X 1 and X 3 with X 0 as given by one-loop dressed perturbation theory. The left plot refers to the inversion using H 1 (q 2 ) = 1, while the right plot are the results for the inversion using the improved gluon-ghost vertex. Smaller values of the regularizing parameter are associated to solutions with larger norms, while larger values of produce X 1 and X 3 with smaller norms. Recall that X 1 has mass dimensions, while X 3 is dimensionless.
equations, is among the smallest values just before the norms of X 1 and X 3 start to grow but without changing the residuum. On the above figure we point out three solutions in the region where takes approximately its optimum value.
Our first comment on Fig. 13 being that both the scalar and vector components of the Dyson-Schwinger equations can be resolved with the ansatz considered, i.e. setting X 0 (p 2 ) to its one-loop dressed perturbative result and getting X 1 (p 2 ) and X 3 (p 2 ) from solving the gap equations, provided we let the norm of X 1 and X 3 be large enough. Of course, for large norms X 1 and X 3 are free to vary over a large range of values and, therefore, solutions with smaller norms are preferred.
From Fig. 13 three typical solution close to the optimal solution, as defined previously, are identified. For the X 0 perturbative solution using the tree level (TL) ghost-gluon vertex, the characteristics of these solutions are have the same m b.m. and Z 2 as the previous ones and ||X 0 − 1|| 2 = 0.45283. The norms of X 0 , X 1 and X 3 are the norms of the vectors that appear in the linear systems. The quality of the solutions can be appreciated on Fig. 14 where we show both the l.h.s. of the scalar and vector components of the gap equation, together with the difference between the l.h.s. and the computed r.h.s. using the X 0 from one-loop perturbation theory and X 1 and X 3 that solve the modified linear system. The relative error both for the scalar and vector components of the Dyson-Schwinger equations are shown on Fig. 15. On the figures we have defined ||∆Sca|| 2 and ||∆Vec|| 2 should be understood as the sum of the squares of the components of the linear systems (6.6) and (6.6), respectively, over the Gauss-Legendre points. As Fig. 15 shows, the relative error on the DSE equations is below 10% for the scalar equation and below 8% for the vector equation. Surprisingly, despite the larger values of ||∆Vec|| 2 relative to ||∆Sca|| 2 , the vector component of the gap equation is better resolved. This is also due to the fact that A(p 2 ) spans a narrower range of values relative to B(p 2 ). We have tried to rescale the linear system by 1/A(p 2 ) for the vector equation and by 1/B(p 2 ) for the scalar equation to try to improve the quality of the solutions, specially at large momentum. However, the numeric solutions of the rescaled linear systems produced X 1 and X 3 that don't seem reasonable and, for example, result in a X 1 at the cutoff that is far away from zero. Further, for the rescaled systems the ∆Sca and ∆Vec are larger than the ones obtained without rescaling the linear system. For all these reasons we disregard the rescaled linear system solutions.
The quark-ghost kernel form factors X 1 and X 3 computed for the various associated to the solutions I (TL) -III (TL) and I (Enh) -III (Enh) can be seen on Fig 16. For X 1 the outcome of resolving the integral equations using either the tree level or the enhanced ghostgluon vertex result on essentially the same function. Further, the various solutions provide essentially the same X 1 (p 2 ), with the exception of III (TL) and III (Enh) that return a suppressed form factor relative to the other solutions. For X 3 the situation is similar, with the form factor associated to the solutions I (TL) and I (Enh) being enhanced at momentum above 2 GeV. Looking at Fig. 15, one can observe that solutions II (TL) and II (Enh) are those with smaller relative errors over the full range of momentum considered. So, from now on we will take these solutions as our best solutions associated to the perturbative X 0 form factor. Note that the scalar equation is solved with a relative error 8% and the vector equation is solved with a relative error 6%.

Solving the Dyson-Schwinger Equations
Let us now discuss the simultaneous computation of X 0 , X 1 and X 3 from the modified linear system of equations that replace the original Dyson-Schwinger integral equations. The procedure to build the linear system as well as the regularization of the corresponding linear system of equations follow the steps described on the previous section.    Figure 16: The quark-ghost kernel form factors X 1 and X 3 computed for the tree level ghost-gluon vertex (left) and the enhanced gluon-ghost vertex (right).

Sol. I (Enh) [GeV] Sol. II (Enh) [GeV] Sol. III (Enh) [GeV]
large linear system as follows that again we refer, as a short name, by B = N X. The upper most component of the large vector X contains the form factor X 0 defined in all the set of Gauss-Legendre points used in the integration over the loop momenta. The remaining components of the large X vector are the form factors X 1 and X 3 defined at the lower first half set of Gauss-Legendre points used in the integration over the momentum. This means that the solution of the linear system (6.8) returns X 0 (p 2 ) for p ∈ [0, Λ] and X 1 (p 2 ) and X 3 (p 2 ) for p ∈ [0, Λ/2]. For X 1 and X 3 and for p > Λ/2 the form factors will be assumed to vanish. In order to fulfil the boundary conditions for X 0 (q 2 ) we write X 0 (q 2 ) = 1 +X 0 (q 2 ) and solve the linear system forX 0 (q 2 ), rebuilding X 0 (q 2 ) at the end. The resulting linear system is then regularized using Tikhonov regularisation scheme and the corresponding N T B = (N T N + )X linear system solved for various . The choice of the optimal regularization parameter will be made following the criteria discussed in Sec. 6.1. We have checked that by interchanging the roles of X 0 , X 1 and X 3 when building the large linear system the solutions are unchanged; see more on that below. The differences only occur for those functions calculated only for q ∈ [0, Λ/2], compared to the version of the linear system were they are computed in the range q ∈ [0, Λ]. In the first case, i.e. for the solutions computed only for q ∈ [0, Λ/2], there appears a discontinuity at q = Λ/2 (recall that the form factors are set to zero for momenta above Λ/2) but for smaller q the form factors of all versions of the linear system are indistinguishable. On Fig. 17  On Fig. 17 we identify four solutions associated to an around its optimal value and whose characteristics are The relative errors for the solutions I -IV of the regularized linear system are show on Fig. 18. In general the solution for the vector component of the equation is satisfactory, with the scalar component of the equation being more demanding and not all of the solutions I -IV resolve the scalar part of the gap equation with a relative error below 10%. Only solutions III and IV resolve the DSE equations with a relative error below 8%. In particular for these solutions the value of ||X 0 − 1|| is of the order of 10 −2 suggesting that the nonperturbative solution prefers having a X 0 1 and, in this sense and for this form factor, are close to the result from perturbation theory discussed on Sec. 6.1. The observed growth of the relative error for p 10 GeV is probably related also to the missing components of X 1 and X 3 which are set to zero for this range of momenta. The form factors X 0 , X 1 and X 3 associated to the solutions I -IV can be seen on Figs. 19 -21. On Fig. 19 besides solutions I -IV we also show the perturbative X 0 (q 2 ) computed using one-loop dressed perturbation theory using both the treel level ghostgluon vertex and its enhanced version. The perturbative solutions and those obtained solving the Dyson-Schwinger equations have rather different structures, with perturbation theory providing larger X 0 (p 2 ) and predicting a relatively large tail. Indeed, the solutions of the regularized linear system recover their tree level value X 0 (p 2 ) = 1 from p 10 GeV onwards, while the perturbative solution only reproduces its tree level value at much larger momentum. The momentum scale associated to the absolute maxima of X 0 (p 2 ) occurs at essentially the same p ≈ 400 MeV, while the perturbative results points to a maximum of X 0 (p 2 ) at momenta slightly above the GeV scale. Qualitatively, the nonperturbative solutions all have the same pattern for this form factor. The exception being Sol. I which clearly overestimates |X 0 | for p 1 GeV. The solutions III and IV are those which resolve the gap equation with the smaller relative errors below 8% -see Fig. 18. The non-perturbative solution of the DSE gives a X 0 (p 2 ) that differs from its tree level value by less than 5%, that are above unit for momenta p 1 Gev. At this momenta scale the form factors reaches values below one, reaching a minimum for p just above 1 GeV, and approaching its tree level value at high momentum from below. The differences between the non-perturbative X 0 and its tree level value for p 10 GeV are rather small.
Our non-perturbative estimations for X 1 (p 2 ) can be viewed on Fig. 20. All the solutions I -IV reproduce the same pattern for this form factors, with a positive maxima around p 400 MeV and with X 1 becoming small for p 1.5 GeV. In particular, for the solutions III and IV, X 1 (p 2 ) is particularly small ( 0.4 GeV) for p 1.5 GeV. One should not forget that the quark-ghost form factor appearing in the quark-ghost kernel is not X 1 but this function times the gluon propagator -see Eq. (4.21). The same applies to X 3 as can be 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Figure 19: X 0 (p 2 ) from inverting the Dyson-Schwinger equations together with its estimation using one-loop dressed perturbation theory by solving exactly Eq. (6.3).  The form factor X 3 (p 2 ) is reported on Fig. 21. It turns out that this function is positive for p 400 MeV and for p 1 GeV, takes negative values in between, takes its maximum value at p 1.5 GeV and then slowly approaches its tree level value from above. Given the way the solutions are computed, on Fig. 21 X 3 (p 2 ) shows a jump at p 10 GeV which corresponds to Λ/2. A similar behaviour can be seen on Fig. 20 for X 1 (p 2 ). However, given that for p 10 GeV X 1 (p 2 ) 0, this sudden jump is not so easily observed.
Finally, on Fig. 22 we provide the various solutions for X 0 (p 2 ), X 1 (p 2 ) and X 3 (p 2 ) after permuting the role of the form factors when writing the extended vector X. For the so-called X 0 X 1 X 3 the extended vector included X 0 over the full set of Gauss-Legendre points with X 1 and X 3 being obtained only in the range q ∈ [0 Λ/2]. For the so-called X 1 X 0 X 3 the extended vector included X 1 over the full set of Gauss-Legendre points with X 0 and X 3 being obtained only in the range q ∈ [0 Λ/2]. For the so-called X 3 X 0 X 1 the extended vector included X 3 over the full set of Gauss-Legendre points with X 0 and X 1 being obtained only in the range q ∈ [0 Λ/2]. We call the readers attention to the stability of the solution of the various linear systems. Furthermore, the comparison of Fig. 16 from Sec.6.1 and Fig. 22 show quite similar X 1 and X 3 suggesting, once again, that X 0 almost does not deviates from its tree level value.

Solving the DSE for
The non-perturbative solutions of the Dyson-Schwinger equations discussed on the previous paragraph suggest that X 0 (p 2 ) 1. Therefore, herein we present the results found when solving the DSE setting X 0 (p 2 ) = 1. The residua of the scalar and vector equations against the norm of the two remaining form factor X 1 and X 3 can be seen on Fig. 23. The characteristics of the solutions highlighted in the figure and associated to an close to its optimal value are  Fig. 24 which shows that solution I resolves the DSE up to p 10 GeV with an error that is smaller than 3% for the scalar equation and error of about 1% for the vector equation. However, the form factor X 3 associated to solution I does not seem to be converged for momenta above 2 GeV. The solution named IV solves the scalar equation with a relative error below 10% and the vector equation with a relative error below 6%.
The form factors X 1 (p 2 ) and X 3 (p 2 ) associated to the solutions I -IV are reported on Fig. 25 and reproduce the same patterns as those of the solutions computed in the previous sections.  Figure 25: The solutions of the DSE for X 1 and X 3 when X 0 = 1.

Full Form Factors and Comparison of Solutions
In Secs. 6.1, 6.2 and 6.3 we have solved the Dyson-Schwinger equations assuming that the quark-ghost kernel form factors are given by Eqs. (4.20) -(4.22). So, besides, X 0 , the full form factors appearing in H and H, see Eqs. (2.16), should be multiplied by the gluon propagator at the proper kinematical configuration. Herein, we aim to compare the various solutions found in previous section and, in this way, provide an estimation of the systematics associated to our ansatz, and also to provide the form of the full functions appearing on H and H. Looking at the relative errors on the Dyson-Schwinger equations and at the convergence of the form factors at higher momenta, the comparison will be done using Sol. II computed using the perturbative X 0 and the tree level ghost-gluon vertex, Sol. III computed when the gap equation is solved for the full set of form factors and Sol. IV when the gap equation is solved for X 0 = 1. Let us start with the X 0 that we have assumed to be only a function of the gluon momenta. The perturbative solutions are compared with the solution obtained inverting the Dyson-Schwinger equations for the full set of form factors used in our ansatz can be seen in Fig. 26. This figure repeats partially Fig. 19 providing a clear view of the solutions. All solutions show a X 0 that essentially is close to its tree level value, i.e. X 0 = 1, while the perturbative solutions have the largest deviation from unit.
The form factor X 1 (p 2 ) can be seen on Fig. 27 for all the solutions. Note that all solutions reproduce essentially the same function of the gluon momentum, with X 1 (p 2 ) being small for p 1.5 GeV and showing a sharp peak at p 400 MeV. X 1 (p 2 ) is positive defined except for a small range of momenta p ∈ [0.75 , 1.4] GeV where it takes small negative values.
The form factor X 3 (p 2 ) can be seen on Fig. 28 for all the solutions. Surprisingly, X 3 (p 2 ) seems to have a relative large tail that appears for all the solutions. Up to momenta p 3 GeV the solutions reproduce essentially the same function. However for p 3 GeV the solution associated to X 0 = 1 is enhanced relative to all the others, with the solutions associated to the one-loop perturbative X 0 being slightly enhanced relative to the nonperturbative solution obtained from inverting the gap equation. X 3 (p 2 ) shows a maxima at p 200 MeV, an absolute maxima at p 1.4 GeV and an absolute minima at p 650 MeV. This form factor is positive defined at infrared momenta p 350 MeV and the high momenta p 900 MeV taking also negative values in p ∈ [0.35 , 0.9] GeV.

Tunning α s
The results for the relative errors on the scalar and vector components of the Dyson-Schwinger equation seen on Figs. 15, 18, 24 show a relative error that for p 10 GeV grow with p and take its maximum value ∼ 10% at the cutoff. This can be viewed in many ways  Figure 28: The quark-ghost kernel form factor X 3 (p 2 ). and one of them being that our choice for the strong coupling constant is not the best one. In our approach we mix quenched lattice results with dynamical simulations and, in order to be able to solve the gap equation for the quark-ghost kernel, the renormalization constant Z 1 , see Eq, (2.6), is set to identity. Although the original integral equation is linear on the form factors X 0 , X 1 and X 3 , the regularized system that is solved introduces an extra parameter that needs to be fixed in the way described above and, therefore, changing the strong coupling constant changes the balance between the regularizating parameter and the various form factors, allowing for adjustments on the solutions. Therefore, the relative errors on the integral equations can be adjusted by changing the strong coupling constant.
In this section, we report on the results of solving the regularized linear system of equations that replace the original equations in the way it is described on Sec. 6. As the figures shows, lowering the value of α s (µ) solves the problem of the increase of the relative error observed in Sec. 6.2. Moreover, of the various solutions considered, for α s (µ) = 0.22 one can observe solutions whose relative error is of the order of ∼ 1% for the scalar equation and ∼ 3-4% for the vector components, the solutions named Sol. I and II in Fig. 30. The relative error associated to the remaining solutions given on Figs. 24, 29, 30 and 31 are larger and, therefore, we take α s (µ) = 0.22 as the optimal value for the strong coupling constant within our approach. The corresponding quarkghost kernel can be seen on Figs. 32, 33 and 34, together with the corresponding solution computed using α s (µ) = 0.295. The solutions for the two values of α s are similar, although those associated to the smaller value of α s achieve higher values. If at momenta p 1 GeV Sol. I takes absolute values that are higher than those of Sol. II, at lower momenta the difference between the two solutions is marginal.

The Quark-Gluon Vertex Form Factors
In the previous section we have computed the quark-ghost kernel form factors X 0 (q 2 ), X 1 (q 2 ), X 3 (q 2 ) that, together with the gluon propagator, define the full form factors as given in Eqs. (4.20), (4.21) and (4.22). Once the full quark-ghost kernel form factors are known, then the longitudinal quark-gluon form factors can be computed using Eqs. (3.7) -(3.10), after performing the rotation to the Euclidean space and identifying the g i (p 2 1 , p 2 2 ) functions to g 0 (p 2 1 , p 2 2 ) = 1 and g 1 (p 2 1 , p 2 2 ) = g 2 (p 2 1 , p 2 2 ) = ∆ For completeness, we write the full expressions for the longitudinal form factors in Euclidean space  Figure 34: The quark-ghost kernel form factor X 3 (p 2 ) computed using α s (µ) = 0.22. Note that by taking into account structures of the quark-ghost kernel other than X 0 the quark-gluon vertex deviates considerable from a Ball-Chiu type and it is now a function both of p, q and of the angle between the quark and gluon momenta. The angular dependence appears associated to the scalar products (pq) and also on the argument of the gluon propagator ∆ (p 2 + k 2 )/2 .
For the calculation of λ 1 -λ 4 we will use Sol. II computed using α s (µ) = 0.22; see Sec. 6.5 for details. We recall the reader that the calculation performed here considers only the longitudinal form factors and that the ansatz for the vertex takes into account the dependence between the angle of the incoming quark momentum and the incoming gluon momentum.
The overall picture of the various form factors when the angle between the incoming quark momentum p and the incoming gluon momentum q is θ = 0 can be seen on Fig. 35. On Fig. 36 the λ 1 to λ 4 are given for a θ = 2π/3. The form factors λ 1 to λ 4 are finite for all p and q and approach asymptotically their perturbative values. Further, for our definition of the operators L µ , see Eqs. (2.12) for their definition in Minkowsky space, the corresponding form factors are essentially positive defined. The exception being λ 4 that takes both positive and negative values and whose maximum absolute value is negative and appears for small p and q. The relative magnitude of the λ i suggest that the quark-gluon vertex is essentially saturated by λ 1 and λ 3 , with λ 2 and λ 4 playing minor roles, i.e. the tensor structures of the longitudinal part of the vertex seem to play a subleading role; see also the discussion for the soft quark limit, defined by a vanishing quark momentum, and the symmetric limit below.
Our result differs significantly from the perturbative estimation of the form factors [5], where all the strength appears associated to λ 1 . For example, for the kinematical configuration defined by p 2 = (p − q) 2 at vanishing p we have λ 1 ≈ 1.1 and λ 2 ≈ 0.12 GeV −2 and λ 3 ≈ 0.18 GeV −1 for a current mass m q = 115 MeV, a renormalization scale µ = 2 GeV and for α s = 0.118. Of course, one should look to the relative values of the various λs and not to their absolute values. For the comparison of the contributions from the various form factors one can use the non-perturbative momentum scale of 1 GeV to build dimensionless quantities. Then, as seen on Figs. 35 and 36 the scales for λ 1 and λ 3 are similar, while the maximum of λ 2 is about 10% relative to the maxima of λ 1 and λ 3 and the maximum for λ 4 is about half of that for λ 2 . The comparison of our results with those reported in [17,22,23] is difficult to perform but in these works λ 1 clearly dominates. On [17] λ 2 reaches at most 16% of the maximum value of λ 1 , while λ 3 seems to have the possibility of taking large values. On [22,23], λ 2 and λ 3 take, at most, a numerical value that is about 23% of the maxima of λ 1 , with λ 4 being essentially negligible. Our solution shows a vertex dominated by λ 1 and λ 3 with these form factors reaching numerical values of the same order of magnitude -see also Fig. 41.  As seen on Figs. 35 and 36 the quark-gluon form factors are significantly enhanced for low values of p and q. The momentum region where one observes the enhancement of the λ 1 to λ 4 happens for p 1 GeV and q 1 GeV, with its maximum values showing up for p ≈ q ≈ Λ QCD -see, also, the discussion below on the angular dependence.
The infrared enhancement of λ 1 to λ 4 with the gluon momentum is a direct consequence of using the Slavnov-Taylor identity (2.14) to rewrite the form factor. Indeed, as can be seen on Eqs. (7.2) -(7.2), all the form factors have, as a global factor, the ghost dressing function F (q 2 ). The ghost dressing function is enhanced in the infrared, see Fig. 4, implying the increase of the λ i as q −→ 0.
The infrared enhanced of the form factors with the quark incoming momentum is more subtle. It is linked to the ansatz considered herein, that relies on the analysis of the soft gluon limit of the Landau gauge lattice data for λ 1 performed in [29]. Indeed, this work identified a dependence of λ 1 on the gluon propagator that was incorporated in our ansatz when the quark-ghost kernel form factors X 1 and X 3 were made proportional to ∆((p 2 + (p − q) 2 /2). This term is crucial to have well behaved kernels in the integral equations, i.e. to ensure that the Dyson-Schwinger equations are finite, and it introduces an additional dependence on the angle between the quark and the gluon momenta. The gluon propagator is a decreasing function of its argument and, therefore, for a given q and angle between the quark and gluon momenta, the terms proportional to X 1 and X 3 increase as p decreases. This explains, in part, the observed enhanced of the quark-gluon form factors.
The quark-gluon form factors are functions of the p, q and of the angle between the two vectors. Their dependence on the angle can be seen on Figs. 37 -40. These figures also provide a clear picture of the maxima of the various form factors as functions of the gluon momenta. For λ 1 and λ 3 the maxima are for q ≈ 300 MeV, while for λ 2 the maximum is at q ≈ 600 MeV. λ 4 seems to be a more complicated function of p, q and θ. Indeed, this later form factor shows various maxima of the same order of magnitude for different p, q and θ values.
All the form factors appear to be monotonous decreasing functions of the angle between the incoming quark and incoming gluon momenta θ. If the pattern of the q dependence of     λ 1 , λ 2 and λ 3 seems to be independent of θ, λ 4 seems to reverse is behaviour relative to the q − axis for θ π/3. Clearly, the maximum values for all the form factors occurs for θ = 0, i.e. the quark-gluon vertex favours the kinematical configurations with small values of p and q and also of the angle between the quark and gluon momentum 1 . It follows that the quark-gluon vertex seems to favour low values of the quark and gluon momenta and a p and a q that are preferably parallel vectors. From the point of view of the momentum dependence, our solution for the quarkgluon vertex is closer to that of the Maris-Tandy model [32] than those computed in [17,22,23]. Recall that the Maris-Tandy model considers a single form factor, that would be equivalent to our λ 1 , and ignores the dependence of the vertex on the quark momentum. In particular, for this model we also checked that the region where the quark-gluon form factors computed here are enhanced occurs essentially within the same range of momenta as for the corresponding form factor of the Maris-Tandy model. Note also that the maxima of the form factors computed in the present work occur for momenta where the kernels appearing in the original equations take there maximum values -see Figs. 6, 9 and 10. It is difficult to measure the relative importance of the contribution of the longitudinal form factors λ 1λ 4 to the quark-gluon vertex. However, an idea of their relative importance can be "measured" looking at particular kinematical configurations. Herein we consider the soft quark limit where the incoming quark momentum vanish and the totally symmetric limit where p 2 = q 2 = k 2 and θ = 2π/3. The corresponding form factors multiplied by appropriated powers of momenta to build dimensionless function can be seen on Fig. 41 (computed using the θ = 2π/3 data). If for the symmetric configuration the dominant form factor seems to be λ 1 , for the soft quark limit that role is played by p λ 3 . Note that the maximum of the later is about 1.3 times larger than the maximum of the former. Curiously, the maxima of λ 1 and pλ 3 occur at exactly the same momentum scale p = 310 MeV. As the figure shows it seems that the quark-gluon vertex is dominated by λ 1 and λ 3 with the tensor structures associated to λ 2 and λ 4 playing a minor role.
Finally, let us consider the soft gluon limit whose λ 1 form factor as recently being computed using lattice simulations [27]. The data was investigated in [29] revealing an important contribution to λ 1 linked with the gluon propagator. Our estimation of λ 1 in the soft gluon limit can be seen on Fig. 42 -see the full curve in black. This curve was (arbitrarely) normalized to reproduce the lattice data at 1 GeV of the β = 5.29 and M π = 295 MeV simulation 2 . Clearly, our ansatz underestimates λ 1 in the infrared region. As discussed in [29], in the soft gluon limit where p is the quark incoming momenta, and in our notation Our solutions has X 1 (0) ≈ 0 GeV and X 3 (0) ≈ 0 and, therefore, it underestimates λ 1 (p 2 ) in the infrared region. Note that herein X 1 and X 3 are assumed to be a function only of the gluon momentum and, due to the integration over the gluon momenta q in the Dyson-Schwinger equations, these form factor are multiplied by q 3 and, therefore, the inversion probably is not able to resolve correctly X 1 (q 2 ) and X 3 (q 2 ) in the deep infrared region. If in the calculation of the soft gluon limit one assumes that X 1 (0) deviates from zero by a small quantity, the agreement with the lattice data is considerable improved both in the infrared and in the ultraviolet. This is represented by the two full curves in colour of Fig. 42 where X 1 (0) is set to a small value. The colour curves suggest a X 1 (0) ∼ 0.5 -0.7 GeV. Further, the agreement in the ultraviolet region can also be improved if X 3 (0) assumes small and positive values; recall that X 3 (q 2 ) approaches zero from the above when q 2 → 0 as can be seen on Fig. 34.

Summary and Conclusions
In this work we investigate the non-perturbative regime of the Landau gauge quark-gluon vertex (QGV), taking into account only its longitudinal components, and relying on lattice results for the quark, gluon and ghost propagators, together with continuum exact relations, namely a Slavnov-Taylor identity and the quark propagator Dyson-Schwinger equation. Furthermore, we incorporate the exact normalisation condition for the quark-ghost kernel form factor X 0 [17]. In addition, we take into account an empirical relation that links the gluon propagator and the soft gluon limit of the form factor λ 1 checked against full QCD lattice simulations [29]. The full set of the quark-ghost kernel tensor structures are taken into account to build an ansatz for the longitudinal quark-gluon vertex that is a function of both the incoming quark p and gluon q momenta, and the angle between p and q.
The quark-ghost kernel requires four scalar form factors X 0 , X 1 , X 2 , X 3 [35]. For the construction of the quark-ghost kernel a perfect symmetry between incoming and outgoing quark momentum is assumed, which simplified the description of the QGV in terms of X 0 , X 1 = X 2 and X 3 . Charge conjugation demands that for the soft gluon limit, defined by q = 0, λ 4 = 0 and our construction implements such constraint. Noteworthy to mention that our ansatz goes beyond the Ball-Chiu type of vertex [6] and includes it as a particular case, when X 1 = X 3 = 0 and X 0 = 1.
The Dyson-Schwinger equations are solved for the quark-gluon vertex that are written in terms of the unknown functions X 0 , X 1 and X 3 . From the point of view of the quark-ghost kernel form factors, these are linear integral equations. The corresponding mathematical problem is ill defined and requires regularization in order to obtain a meaningful solution. The original integral equations for the scalar and vector components of the quark gap equation are transformed into a set of linear system using Gauss-Legendre quadratures to perform the integrations and after doing the angular integration. In our approach we rely on the Tikhonov linear regularization that is equivalent to minimize ||B − N X|| 2 + ||X|| 2 . The solutions are found numerically after writing the regularized linear system in its normal form. The small parameter is set by looking at the balance between the associated error on the Dyson-Schwinger equations, i.e. the difference between the l.h.s and the r.h.s. ||B − N X|| 2 , and the norm of the corresponding quark-ghost form factors, i.e. ||X|| 2 , for each solution of the regularized linear system.
The resulting quark-gluon vertex form factors λ 1λ 4 show a strong enhancement in the infrared region and deviate significantly from their tree level results for quark and gluon momenta below ∼ 2 GeV. At high momentum the form factors approach their perturbative values. In what concerns the gluon momentum, the observed infrared enhancement for the QGV form factors can be traced back to the multiplicative contribution of the ghost dressing function introduced through the Slavnov-Taylor identity. Recall that the gluon dressing function peaks at q = 0 and, therefore, favours that the incoming and outgoing quark momentum to be parallel. On the other hand, the infrared enhancement associated to the quark momentum is linked to the gluon dependence that was observed on the analysis of the soft gluon limit of the QGV and, clearly, favours small quark momentum p ∼ 0 and also p parallel to q; see Eqs. (4.21), (4.22) and (7.2) -(7.5) and, in particular, the argument appearing on the gluon propagator term. The maxima of the computed form factors are essentially at the maxima of X 0 , X 1 and X 3 and they appear for momenta of p, q ∼ Λ QCD , which again seems to set the appropriate non-perturbative momentum scale. Recall that the momentum scale comes from the use of lattice data for the propagators. Further, we find that the quark-gluon vertex is dominated by the form factors associated to the tree level vertex γ µ and to the scalar 2 p µ + q µ , with the higher rank tensor structures giving small contributions. Overall, our findings are in qualitative agreement with previous works both with phenomenological approaches, as in the case of the Maris-Tandy vertex [32], and those based on first principles ab initio continuum methods, see e.g. [23] and references therein.
The high momentum behaviour of the quark-gluon vertex form factors reproduces their perturbative values. However, the matching between the computed form factors and their perturbative tail is not yet implement. In addition, we verified that for the soft gluon limit, λ 1 is not able to reproduce quantitatively the lattice data from full QCD simulations, apart the qualitative momentum behaviour. This can be traced back to the poor resolution of the kernel in the deep infrared region, due to the q 3 factor coming from the momentum integration. As we have verified, a small tuning of X 1 and X 3 at q = 0 is enough to reproduce the soft gluon limit lattice data within the present framework. This two challenging problems, together with inclusion of the transverse part of the vertex, call for an improvement of the approach devised herein and are to be tackled in a future work. Despite of that, we expect that the present results can help understanding the non-perturbative dynamics of quarks and gluons in the infrared region and that can motivate further applications to the study of hadron phenomenology based on quantum field theoretical approaches as those using Bethe-Salpeter and/or Faddeev equations -see e.g. [45] and references therein.

A 4D Spherical Coordinates and integration over momentum
In 4D the spherical coordinates are related to the cartesian coordinates as follows where Λ stands for the cutoff introduced to regulate the theory.

B Comparing Propagator Fits with Previous Works
For completeness and in order to allow for a better comparison between the of the current work with those reported in [28], we provide the fits used in both works with the gluon propagator, the ghost propagator and the quark wave function curves renormalized at µ = 4.3 GeV within the MOM scheme. On Fig. 43 the curves referred to as JHEP are those of [28], while those designated as NEW are the curves mentioned in Secs. 5.1 and 5.2. As the figure shows, there are differences between the two sets of curves, not only at the infrared region but also on the running at high momentum.