The Quark-Gluon Vertex and the QCD Infrared Dynamics
O. Oliveira W. de Paula T. Frederico J. P. B. C de Melo








e1e-mail: orlando@uc.pt
e2e-mail: wayne@ita.br
e3e-mail: tobias@ita.br
e4e-mail: joao.mello@cruzeirodosul.edu.br

CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
          Dep. de Física, Instituto Tecnológico da Aeronáutica, Centro Técnico Aeroespacial, 12.228-900 São José dos Campos, São Paulo, Brazil
          Laboratorio de Física Teórica e Computacional, Universidade Cruzeiro do Sul, São Paulo, 01506-000 SP, Brazil






The Quark-Gluon Vertex and the QCD Infrared Dynamics
    Orlando Oliveirae1,addr1,addr2
        Wayne de Paulae2,addr2
        Tobias  Fredericoe3,addr2
        J. P. B. C de Meloe4,addr3

    Received: date / Accepted: date
=======================================================================================================================================




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 γ_μ and  to the operator  2   p_μ + q_μ. The higher rank tensor structures provide small contributions to the vertex.

























§ INTRODUCTION


The interaction of quarks and gluons is described by Quantum Chromodynamics <cit.>, 
a renormalisable gauge theory associated to the color gauge group SU(3). 
Of its correlation functions, the quark-gluon vertex has a fundamental role in hadron phenomenology, in the understanding of chiral symmetry breaking 
mechanism and the realisation of confinement. 
Despite its relevance for strong interactions, our knowledge of the quark-gluon 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 <cit.> but only for some 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 what the authors called its asymptotic limit.
In particular, the vertex asymptotic limit was used to investigate  ansätze that can be found in the 
literature <cit.> with the aim to test their 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 <cit.>. 
Typically, the computation is performed after writing the vertex in terms of other QCD vertices and propagators and taking into account its perturbative tail. 
Most of the computations  include only a fraction of the twelve form factors required
to fully describe the quark-gluon vertex. 
In <cit.> the authors look for a first principle determination of the vertex by solving the theory at the non-perturbative level,
gathering information on the vertex from QCD symmetries  and relying on one-loop dressed perturbation theory.
The vertex has also been investigated perturbatively within massive QCD, i.e. using the Curci-Ferrari model <cit.>,
and all its (perturbative) tensor structures form factors accessed for some kinematical configurations.

Lattice simulations, both for quenched <cit.> and full QCD <cit.>, 
were also used to investigate the quark-gluon vertex. Again, only a limited set of kinematical configurations were accessed and, in particular, its the soft gluon limit, 
defined by a vanishing gluon momenta, was mostly explored. 
For full QCD so far only a single form factor, that associated with the tree level tensor structure, was measured on lattice simulation in the soft gluon limit.

One can also find in the literature attempts to combine continuum non-per­tur­ba­ti­ve QCD equations with  lattice simulations
to study the quark-gluon vertex. Indeed, in <cit.> a generalised Ball-Chiu vertex was used in the quark gap equation, together with lattice
results for the quark, gluon and ghost propagators to investigate the quark-gluon vertex. In <cit.>, the full QCD lattice data for
λ_1 was studied 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 equations requires assuming some type of functional dependence for 
various propagator 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. This approach that combines lattice information and continuum equations
relies strongly on the effort to access high precision lattice simulations.

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 
renormalisation <cit.>.
A popular and quite successful model was set in <cit.>, named the Maris-Tandy model,
that assumes a bare quark-gluon vertex and introduces an effective gluon propagator that is strongly enhanced at infrared scales and recovers
the one-loop behaviour at higher momentum. This model simplifies considerably the momentum dependence of the combined 
effective gluon propagator and quark-gluon vertex and assumes that the dominant momentum dependence is associated only with the gluon momentum.
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 effective gluon propagator.
Although the Maris-Tandy model is quite successful for phenomenology, it is not able to describe the full set of hadronic properties 
and fails to explain the mass splittings of the ρ and a_1 parity partners, underestimates 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, <cit.> and references therein. 
Several authors have tried to improve the Maris-Tandy model either by studying its dependence on the various parameters, see e.g. <cit.>, or
by changing its functional dependence at low momenta, see e.g. <cit.>, to achieve either a better description of Nature or good agreement with
the results from lattice simulations.

The goal of the present work is to explore further the quark-gluon vertex in the non-perturbative regime from first principles calculations 
combining continuum methods with results coming from lattice simulations.
Our approach follows the spirit of the calculation performed in <cit.> that solves the quark Dyson-Schwinger equation for the vertex.
In <cit.> the quark-gluon vertex was described as a generalised Ball-Chiu vertex and single unknown form factor,
function only of the gluon momentum, was considered. The current work goes beyond this approximation and includes the full set of longitudinal form factors
that appear in the Ball-Chiu vertex. 
In this work we disregard any contribution due to the transverse form factors and 
consider the Landau gauge quark-gluon vertex, to profit from the recent high quality lattice data for the quark, the gluon and the ghost propagators.
Moreover, our 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, for the infrared vertex, associated with the gluon propagator <cit.>.
As in other studies, we rely on a Slavnov-Taylor identity to write the vertex longitudinal form factors as a function of the quark wave function, 
the running quark mass, the quark-ghost kernel form factors and the ghost propagator. 
The normalisation of the quark-ghost kernel form factors X_0 (see below for definitions) derived in <cit.> for the soft gluon limit is also
taken into account when solving the quark gap equation for the vertex.
The normalisation of X_0 played 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 γ_μ, performed in <cit.> that identified
an important contribution for λ_1 coming from the gluon propagator.

Our solution for the quark-gluon vertex returns a X_0 that deviates only slightly from the normalisation condition referred above.
However, the longitudinal form factors describing the quark-gluon vertex 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 can be traced back to
ghost contribution introduced by the Slavnov-Taylor identity and the gluon dependence of the ansatz. 
At high momentum the form factors seem to approach their perturbative values. 
The matching with the perturbative tail is not perfect and this result can be understood partially due to the regularisation for the mathematical problem,
i.e. the Tikhonov regularisation, and partially to the parametrisation of the vertex.
Indeed, by calling the gluon propagator to describe the various form factors, the inversion of the Dyson-Schwinger equations is
quite sensitive to the low momentum scales, where the gluon propagator is much larger, and less sensitive to the ultraviolet regime. 
In order to overcome this problem, we considered a relatively large cutoff in the inversion and in this way add information on the perturbative tail.

The computed quark-gluon vertex is 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.
The enhancement occurs preferably at momenta of ∼ Λ_QCD.
Furthermore, the vertex is enhanced when all momenta entering the vertex (see Fig. 1) tends to be parallel  in pairs,  
solving in this way the compromise  that the momenta are restricted to a region around Λ_QCD.
Within our solution for the quark-gluon vertex, the dominant form factors are associated with the tree level vertex γ_μ and the operator 2   p_μ + q_μ.
The higher rank tensor structures give sub-leading contributions to the vertex.

In the current work, the vertex is written using the Ball-Chiu construction. It is known that the Ball-Chiu vertex has kinematical
singularities for the Landau gauge <cit.> that are associated to the transverse form factors (see definitions below).
These singularities can be avoided by considering a different tensor basis for the full vertex as described, for example, in <cit.>. 
However, the singularities are not associated to the longitudinal form factors and our calculation only takes into account this class of form factors.

The paper is organised as follows. In Sec. <ref> 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 and the quark-ghost kernel. The parametrisation
of the quark-ghost kernel is also discussed. In Sec. <ref> the scalar and vector components of the DSE in Minkowski space are given,
together with the corresponding kernels.
In Sec. <ref> the DSE are rewritten in Euclidean space, introduce the vertex  ansatz and 
perform a scaling analysis of the integral equations. 
In Sec. <ref> we give the details of the lattice data used in the current work for the various propagators
and on the functions that parametrise the lattice data. The kernels for the Euclidean space DSE are discussed in Sec. <ref>,
together with the solutions for the vertex of the gap equation. The quark-gluon vertex form factors are reported in Sec. <ref> for several
kinematical configurations.
Finally, on Sec. <ref> 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. In this first part of this work, the equations discussed 
are written in Minkowski space with the diagonal metric g = ( 1,   -1,   -1,   -1).
Let us follow the notation of <cit.> for the quark-gluon vertex represented in  Fig. <ref>
that considers all momenta are incoming and, therefore, verify

    p_1+p_2+p_3 = 0  .

The one-particle irreducible Green's function associated to the vertex reads

    Γ^a_μ (p_1, p_2, p_3) = g   t^a  Γ_μ (p_1, p_2, p_3)  ,

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

    S(p)     =    i/A(p^2)   p -  B(p^2) 
      = i  A(p^2)   p + B(p^2) / A^2(p^2)    p^2 - B^2(p^2)
        =    
      i   Z(p^2)    p + M(p^2) / p^2 - M^2(p^2)  ,

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 renormalisation group invariant running
quark mass. 

The Dyson-Schwinger equation for the quark propagator, also named the quark gap equation, is represented in Fig. <ref>
and can be written as

    S^-1(p)  =   -i   Z_2 (   p - m^bm) + Σ (p^2)  ,

where Z_2 is the quark renormalisation constant, m^bm the bare current quark mass and the quark self-energy is given by

    Σ (p^2)     =     Z_1  ∫d^4 q(2π)^4 ×
          × D_μν^a b(q)    (  i  g  t^b γ_ν )   S(p-q)  Γ^a_μ (-p,   p-q,   q) ,

where Z_1 is a combination of several renormalisation constants, D_μν^a b(q) is the gluon propagator that, in the Landau gauge,
is given by

    D^ab_μν (q) =  -i δ^ab( g_μν - q_μ q_ν/q^2) D ( q^2 )  ;
   
below both D^ab_μν (q) and 
D (q^2) will be referred to as the gluon propagator.




A key ingredient in gap equation (<ref>) is the quark-gluon vertex. Indeed, it is only after knowing Γ^a_μ or, equivalently Γ_μ, that 
Z(p^2) and M(p^2) can be computed.
The Lorentz structure of the quark-gluon vertex Γ_μ can be decomposed into 
longitudinal Γ^(L) and transverse Γ^(T) components relative to the gluon momenta, i.e. one writes

    Γ_μ(p_1,   p_2,   p_3)  = Γ^(L)_μ(p_1,   p_2,   p_3) + Γ^(T)_μ(p_1,   p_2,   p_3),

where, by definition,  

    p_3^μ Γ^(T)_μ(p_1,   p_2,   p_3)= 0.

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  form factors and
for the Ball and Chiu basis <cit.> it reads

    Γ^L_μ (p_1, p_2, p_3)     =     -i ∑^4_i=1λ_i (p_1, p_2, p_3)   L^(i)_μ (p_1 , p_2)  
    
       Γ^T_μ (p_1, p_2, p_3)      =     -i ∑^8_i=1τ_i (p_1, p_2, p_3)   T^(i)_μ (p_1 , p_2)  .

The operators associated to the  longitudinal vertex are

    L^(1)_μ (p_1 , p_2)     =     γ_μ , 
    
     L^(2)_μ (p_1 , p_2)     =      (  p_1 -  p_2 ) ( p_1 - p_2 )_μ , 
    
     L^(3)_μ (p_1 , p_2)     =      ( p_1 - p_2  )_μ 𝕀_D  ,  
    
     L^(4)_μ (p_1 , p_2)     =    σ_μν( p_1 - p_2 )^ν ,

while those associated to the transverse part of the vertex read

    T^(1)_μ (p_1 , p_2)      =     [ p_1μ(p_2 · p_3 ) -  p_2μ(p_1 · p_3 )  ] 𝕀_D   , 
    
    
     T^(2)_μ (p_1 , p_2)      =        - T^(1)_μ (p_1 , p_2) (  p_1 -  p_2 )  ,  
    
    
     T^(3)_μ (p_1 , p_2)      =      p^2_3  γ_μ - p_3μ  p_3  ,


    T^(4)_μ (p_1 , p_2)      =         T^(1)_μ  (p_1 , p_2)  σ_αβ  p^α_1 p^β_2  , 
    
    
     T^(5)_μ (p_1 , p_2)      =    σ_μν  p^ν_3  ,  
    
    
     T^(6)_μ (p_1 , p_2)      =    γ_μ( p^2_1 - p^2_2 ) + ( p_1 - p_2 )_μ  p_3  , 
    
    
     T^(7)_μ (p_1 , p_2)     =      -1 2( p^2_1 - p^2_2 ) 
                          [ γ_μ(  p_1 -  p_2 ) - ( p_1 - p_2 )_μ𝕀_D ]  
           -
                          ( p_1 - p_2 )_μσ_αβ  p^α_1   p^β_2  , 
    
    
     T^(8)_μ (p_1 , p_2)     =       - γ_μ σ_αβ  p^α_1   p^β_2 + ( p_1μ   p_2 -  p_2μ  p_1)  ,

















where σ_μν = 1 2 [γ_μ,γ_ν]. 





 §.§ 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. <cit.> and references therein. 
On the  other hand,  gauge symmetry implies that the Green functions also satisfy the Slavnov-Taylor identities (STI)  <cit.>.
These 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

    p^μ_3  Γ_μ ( p_1, p_2, p_3 )     =     F(p^2_3)   [ S^-1 ( -p_1)   H( p_1, p_2, p_3 )  
          
        -   H (p_2, p_1, p_3)   S^-1(p_2)  ]  ,

where the ghost-dressing function F(q^2) is related to the ghost two-point correlation function as

    D^a b (q^2) = - δ^ab  D_gh (q^2) = - δ^ab F(q^2)/q^2
 
and H and H are associated to the quark-ghost kernel.
As discussed in <cit.>, these functions can be parametrised in terms of four form factors as

    H( p_1, p_2, p_3 )     =      X_0 𝕀_D +  X_1   p_1 +  X_2   p_2 + X_3  σ_αβ p^α_1 p^β_2  ,  
    H( p_2, p_1, p_3 )     =    X_0 𝕀_D -  X_2  p_1 - X_1  p_2 +
         X_3  σ_αβ p^α_1 p^β_2   ,

where X_i ≡ X_i ( p_1, p_2, p_3 ) and X_i ≡ X_i( p_2, p_1, p_3 ). 

The STI given in Eq. (<ref>) can be solved with respect to the vertex <cit.> to write the longitudinal form factors λ_i  
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

    λ_1 (p_1, p_2, p_3)  =  F(p^2_3)/2{
           A(p^2_1) [ X_0 + ( p^2_1 - p_1 · p_2 ) X_3 ] 
          
                  +     A(p^2_2) [ X_0 + ( p^2_2 - p_1 · p_2 ) X_3 ]  
          
                  +   B(p^2_1) [X_1 + X_2 ]     
           +   B(p^2_2) [ X_1 + X_2 ]    } ,
                  
          λ_2 (p_1, p_2, p_3)  = 
                   F(p^2_3)/2( p^2_2 - p^2_1 ) {
           A(p^2_1)[ ( p^2_1 + p_1· p_2 ) X_3  - X_0 ] 
                   
          
                   +     A(p^2_2) [ X_0 - ( p^2_2 + p_1· p_2 ) X_3 ]  
          
                   +    B(p^2_1) [X_2 - X_1 ]  
           
                   +    B(p^2_2) [ X_1 - X_2 ]    } ,
                    
          λ_3 (p_1, p_2, p_3)  = 
                       F(p^2_3)/p^2_1 - p^2_2{
           A(p^2_1) [ p^2_1   X_1 + p_1 · p_2  X_2] 
          
                       -     A(p^2_2) [ p^2_2 X_1 + p_1 · p_2  X_2 ]   
          
                          +     B(p^2_1)   X_0 
          
                          -    B(p^2_2)  X_0    } ,
                    
          λ_4 (p_1, p_2, p_3)  =  - F(p^2_3)/2{ A(p^2_1)   X_2 -  A(p^2_2)  X_2 
          
                           +    B(p^2_1)   X_3  - B(p^2_2)  X_3   } .

A nice feature of the above solution for the various form factors λ_i, that can be checked 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 when modelling the vertex.





§ DECOMPOSING THE DYSON-SCHWINGER EQUATION INTO ITS SCALAR AND VECTOR COMPONENTS 


The Dyson-Schwinger equation for the quark propagator is written in (<ref>), with the quark self-energy being given by (<ref>).
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 (<ref>) which, after some algebra, reduces to

    i B(p^2)    =   i Z_2  m^bm
           +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2) /[ A(k^2) ]^2    k^2 - [ B(k^2) ]^2{
           
              2   h(p,q)  ( 2  [ B(k^2) ]  λ_2 -   [ A(k^2) ]  (  λ_3 + λ_4 ) ) 
           +  3 (  [ B(k^2) ]  λ_1 + [ A(k^2) ]  ( 2 p^2 + q^2 - 3 (p ·  q) ) λ_4) 
      } ,

after insertion of the  vertex decomposition (<ref>), taking into account only its longitudinal part,  where k = p -q,

    h(p,q) =  p^2    q^2   -    (p· q)^2 / q^2  ,

 λ_i ≡λ_i (-p ,   p-q,   q) and C_F = 4/3 is the Casimir invariant associated to the SU(3) fundamental representation.

The vector component of (<ref>) is obtained after multiplication by p and then taking the trace of the resulting equation to
arrive on

    -i p^2 A(p^2)   =   -i Z_2  p^2 
           
      +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2) /[ A(k^2) ]^2    k^2 - [ B(k^2) ]^2×
          ×{ 2   h(p,q)  ( [ A(k^2) ]   [ λ_1 +  λ_2  (2p^2 + q^2 - 2 p· q  ) ]  
          
                                                           + [ B(k^2) ] (λ_4 - λ_3) ) 
          
     +   3  (  [ A(k^2) ]  λ_1  ( ( p· q ) - p^2  )  
          
     + [ B(k^2) ]  λ_4   ( (p · q) - 2 p^2  ) ) } .


The two equations (<ref>) and (<ref>) can be simplified further by modelling the quark-gluon vertex.
For example, in <cit.> the vertex was parametrised using the solution of the Slavnov-Taylor identity 
(<ref>)–(<ref>) 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 quark-ghost kernel but X_0, 
assumes that at the non-perturbative level the hierarchy of the form factors follows its relative importance  observed in the high momentum regime. 
Furthermore, in order to compute a solution of the Dyson-Schwinger equations it was introduced a further restriction on X_0, that it depends 
only on the incoming gluon momenta, i.e. that X_0 = X_0 (q^2). 

In order to solve the Dyson-Schwinger equations for the vertex, it will be assumed that the form factors associated to the quark-ghost kernel, see Eq. (<ref>),
factorize as

    X_i (p^2_1, p^2_2, p^2_3) =  g_i(p^2_1, p^2_2)   Y_i (p^2_3)  ,
   
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 description of the quark-gluon vertex <cit.> and simplifies considerably the analysis of the solutions of the equations to be solved.
In <cit.> it was proved that, to all-orders,
 
    X_0 (p,    -p,   0 ) = 1    X_1 (p,    -p,   0 ) = X_2(p,    -p,   0 )

that, for the ansatz (<ref>) implies that

    g(p^2,p^2)   Y_0(0)         =  1  , 
    
    
    
       g_1(p^2, p^2)   Y_1(0)         =  g_2(p^2, p^2)   Y_2(0)  .

A solution that complies with the second relation given in Eq. (<ref>) is to assume 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.
If one takes into account these relations into the ansatz for the quark-gluon vertex, then the solutions of the Slavnov-Taylor identities 
(<ref>)–(<ref>) become

    λ_1 (-p,   p-q,   q )  =  F(q^2)/2{
          [ A(p^2)  + A( k^2 ) ]   g_0(p^2, k^2)   Y_0(q^2) 
           +   2  [ B(p^2)  + B( k^2 ) ]   g_1(p^2, k^2)   Y_1(q^2)  
           +  [ A(p^2) ( 2p^2 - (pq) ) 
           + A(k^2) ( 2 p^2 + q^2 - 3 (pq) )   ]  
                              ×
          ×  g_3(p^2, k^2)   Y_3(q^2)  }  ,




    λ_2 (-p,   p-q,   q ) =  F(q^2)/2( q^2 - 2 (p · q))×
          ×{Δ A      g_0(p^2, k^2)   Y_0(q^2)  
          
                                             + [ A(p^2)   (pq)  - A(k^2) ( q^2 -  (pq) ) ] ×
          ×   g_3 (p^2, k^2)   Y_3(q^2) }  ,




    λ_3 (-p,   p-q,   q ) =  F(q^2)/q^2 - 2 (p · q)×
          ×{Δ B     g_0 (p^2, k^2)    Y_0(q^2) 
           +  [ A( k^2 )  ( q^2 -   (p q) ) - A(p^2)   ( p q )   ] ×
          ×   g_1 (p^2, k^2)   Y_1(q^2) }  ,


    λ_4 (-p,   p-q,   q )  =  F(q^2)/2×
          ×{Δ A    g_1 (p^2, k^2)   Y_1(q^2)  + Δ B    g_3 (p^2, k^2)   Y_3(q^2) }  ,

where

    Δ A = A(k^2) - A(p^2)   , Δ B = B(k^2) - B(p^2)   ,

and k = p - q.
The scalar component of the Dyson-Schwinger equations is now

    i B(p^2)       =     i   Z_2  m^bm
           +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2)   F(q^2)/[ A(k^2) ]^2    k^2 - [ B(k^2) ]^2×
          ×{ g_0(p^2, k^2)   Y_0 (q^2)  𝒦^(0)_B (p,q) 
           +   g_1(p^2, k^2)   Y_1(q^2)  𝒦^(1)_B (p,q) 
          
                                   +    g_3(p^2, k^2)   Y_3(q^2)  𝒦^(3)_B (p,q)  }

where the kernels are defined as

    𝒦^(0)_B (p,q)       =    2   h(p,q)/q^2 - 2(pq) ( B(k^2)  Δ A   -    A(k^2)   Δ B )  
           +  3/2  B(k^2)   [ A(k^2) + A(p^2) ]   ,
   



    𝒦^(1)_B (p,q)    =  
           =   h(p,q)    A(k^2)  [  2  A(p^2) (p q) - A(k^2) (q^2 -   (p q) ) /q^2 - 2 (p q)  - Δ A ]
      
           +   3   B(k^2)  (  B(p^2) + B(k^2)  ) 
            +  3/2  A(k^2)  Δ A  ( 2 p^2 + q^2 - 3(p q) )   ,
  



    𝒦^(3)_B (p,q)    =  
          
     h(p,q)  { 2   B(k^2)    A(p^2) (p q) - A(k^2) (q^2 -   (p q) ) /q^2 - 2 (p q) 
           - A(k^2)  Δ B }
           +  3/2  B(k^2)  [ A(p^2) ( 2p^2 - (pq) ) 
           +   A(k^2) ( 2p^2 + q^2 - 3 (pq) ) ] 
           +  3/2  A(k^2)  Δ B  ( 2 p^2 + q^2 - 3(p q) )    .

Similarly, the vector component of the Dyson-Schwinger equations reduces to

    -i p^2 A(p^2)   =   -i   Z_2  p^2 
           +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2)  F(q^2) /[ A(k^2) ]^2    k^2 - [ B(k^2) ]^2×
          ×{ g_0(p^2, k^2)   Y_0(q^2)  𝒦^(0)_A (p,q) 
           +   g_1(p^2, k^2)   Y_1 (q^2)  𝒦^(1)_A (p,q)  
          
           +    g_3(p^2, k^2)   Y_3 (q^2)  𝒦^(3)_A (p,q) }

with the kernels given by

    𝒦^(0)_A (p,q)  =  
           
       =  h(p,q)  { A(k^2)  [ A(p^2) + A(k^2)  + 2p^2 + q^2 - 2(p q)/q^2 - 2(p q) Δ A ] 
           
                                                    - 2 B(k^2)   Δ B/q^2 - 2 (p  q)}
           - 3/2   A(k^2)    ( A(p^2) + A(k^2) ) (p^2 - ( p q ) )   ,




    𝒦^(1)_A (p,q)  =  
           =  h(p,q) { 2    A(k^2)  ( B(p^2) + B(k^2)  )  
           +   B(k^2) [ Δ A + 2 A(p^2) (p  q)  - A(k^2) (q^2 -  (p q)) /q^2 - 2 (p q)] }
            -   3   A(k^2)  ( B(p^2) + B(k^2)  )    (p^2 - (p  q) )    
            +  3/2  B(k^2)   Δ A  ( (p q) - 2   p^2 )   ,




    𝒦^(3)_A (p,q)  =  
            h(p,q) {  A(k^2)  [ A(p^2) (2p^2 - (pq) )  
           + A(k^2) ( 2p^2 + q^2 - 3 (pq) ) 
             + 2p^2 + q^2 - 2(pq)/q^2 - 2(pq)×
          ×( A(p^2) (pq) - A(k^2) (q^2 - (pq)) ) ] 
           +   B(k^2)    Δ B }
           +  3/2( (pq) - p^2 )   A(k^2) ×
          ×( A(p^2) ( 2p^2 - (pq) ) + A(k^2) ( 2 p^2 + q^2 - 3(pq) ) )  
           +  3/2( (pq) - 2 p^2 )   B(k^2)  Δ B   .






§ THE DYSON-SCHWINGER EQUATIONS IN EUCLIDEAN SPACE 


As already stated, our goal 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 the propagators we will rely on lattice inputs that provide first principles 
non-perturbative 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 

    p^2   → -p_E^2   (p · q)   → - (p_E · q_E)
    A(p^2)   → A_E (-p_E^2)   B(p^2)   → B_E (-p_E^2)
    ∫_q   → i ∫_q_E   Δ(p^2)   → - Δ_E(-p_E^2)

on Eqs. (<ref>)  and (<ref>). For completeness, we provide now all expressions in Euclidean space.

The scalar component of the Dyson-Schwinger equations reads

    B(p^2)   =   Z_2  m^bm
           +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2)   F(q^2)/[ A(k^2) ]^2    k^2 + [ B(k^2) ]^2×
          ×{ g_0(p^2, k^2)   Y_0 (q^2)  𝒦^(0)_B (p,q) 
           +   g_1(p^2, k^2)   Y_1(q^2)  𝒦^(1)_B (p,q) 
          
                                   +   g_3(p^2, k^2)   Y_3(q^2)  𝒦^(3)_B (p,q)  }

while its vector component is given by

    p^2 A(p^2)   =   Z_2  p^2 
           +   C_F Z_1 g^2  ∫d^4 q/(2  π)^4 Δ(q^2)  F(q^2) /[ A(k^2) ]^2    k^2 + [ B(k^2) ]^2×
          ×{ g_0(p^2, k^2)   Y_0 (q^2)  𝒦^(0)_A (p,q) 
           +   g_1(p^2, k^2)   Y_1(q^2)  𝒦^(1)_A (p,q) 
          
                                   +   g_3(p^2, k^2)   Y_3(q^2)  𝒦^(3)_A (p,q)  } .

The kernels appearing in Eqs. (<ref>) and (<ref>) are

    𝒦^(0)_B (p,q)       =    2   h(p,q)/q^2 - 2(pq) [ B(k^2)  Δ A   -    A(k^2)   Δ B ]  
           +  3/2  B(k^2)   [ A(k^2) + A(p^2) ]    ,
  



    𝒦^(1)_B (p,q)    =  
           = 
     h(p,q)    A(k^2)   [  Δ A - 2  A(p^2) (p q) - A(k^2) (q^2 -   (p q) ) /q^2 - 2 (p q) ]
      
           +   3   B(k^2)  [  B(p^2) + B(k^2)  ] 
           +  3/2  A(k^2)  Δ A  ( 3(p q) - 2 p^2 - q^2  )    ,
 



    𝒦^(3)_B (p,q)    = 
           = 
     h(p,q)  { A(k^2)   Δ B 
           - 2   B(k^2)    A(p^2) (p q) - A(k^2) (q^2 -   (p q) ) /q^2 - 2 (p q) }
           +  3/2  B(k^2)  [ A(p^2) (  (pq) - 2p^2  )  
           + A(k^2) (  3 (pq) - 2p^2 - q^2  ) ] 
           +  3/2  A(k^2)  Δ B  ( 3(p q)  - 2 p^2 - q^2 )    .




    𝒦^(0)_A (p,q)  = 
           
        h(p,q)  { -  A(k^2)  [ A(p^2) + A(k^2)  + 2p^2 + q^2 - 2(p q)/q^2 - 2(p q) Δ A ] 
          
                                                    - 2 B(k^2)   Δ B/q^2 - 2 (p  q)}
                                +  3/2   A(k^2)    [ A(p^2) + A(k^2) ] (p^2 - ( p q ) )   ,




    𝒦^(1)_A (p,q)  = 
            h(p,q) { -   2    A(k^2)  ( B(p^2) + B(k^2)  ) 
          
      -   B(k^2) [ Δ A + 2 A(p^2) (p  q)  - A(k^2) (q^2 -  (p q)) /q^2 - 2 (p q)] }
              +   3   A(k^2)  [ B(p^2) + B(k^2)  ]    (p^2 - (p  q) )   
            +  3/2  B(k^2)   Δ A  ( 2   p^2 - (p q)  )   ,




    𝒦^(3)_A (p,q)  =  
           =  h(p,q) {  A(k^2)  [ A(p^2) (2p^2 - (pq) ) 
           +   A(k^2) ( 2p^2 + q^2 - 3 (pq) ) 
             +  2p^2 + q^2 - 2(pq)/q^2 - 2(pq)×
          ( A(p^2) (pq) - A(k^2) (q^2 - (pq)) ) ] 
           -   B(k^2)  Δ B }
           +  3/2( (pq) - p^2 )   A(k^2)  ×
          ×[ A(p^2) ( 2p^2 - (pq) )  
           +   A(k^2) ( 2 p^2 + q^2 - 3(pq) ) ]  
           +  3/2(  2 p^2 - (pq)  )   B(k^2)  Δ B   .


The computation of any solution of the above equations, using lattice inputs for the propagators, requires the use of the renormalised 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)   Y_i (q^2)
and, in particular,  its possible behaviour in the limits where q → 0 and p → + ∞.

Let us start by investigating the ultraviolet limit of the integrand functions appearing 
in Eqs. (<ref>) and (<ref>). In the large q limit it follows that

    q^3   D (q^2)   F(q^2) /A^2(k^2)   k^2 + B^2(k^2)⟶1/q ,

up to logarithmic corrections associated to the various propagators. In this limit, the integrand function appearing on the scalar equation
(<ref>) read

    g_0 (p^2, (p-q)^2 )   Y_0 (q^2)  1/q {3/2 B(q^2) [ A(q^2) + A(p^2) ] }
             +   g_1 (p^2, (p-q)^2 )   Y_1 (q^2)   q  {3/2 A(q^2) [ A(p^2) - A(q^2) ] }
            +   g_3 (p^2, (p-q)^2 )   Y_3 (q^2)    q  ×
          ×{3/2[ A(q^2) B(p^2) - 2   A(q^2) B(q^2)  ] } .

The requirement of having a finite integral demands that at large q

    g_1 (p^2, (p-q)^2 )   Y_1 (q^2)   ≈1/q^2
              g_3 (p^2, (p-q)^2 )   Y_3 (q^2) ≈1/q^2
 
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 
introduced by the renormalisation group analysis are, for large momenta, of type
( log (q^2/ Λ^2) )^γ, with γ standing for 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 (<ref>) times the log correction, the integration function at high momenta becomes

    1/q^2 [ log( q^2/Λ^2_QCD) ]^ -  137/116 ,

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. <ref>, where one can observe the effect due to the log corrections that
suppress further the integrand function at high momenta.

In what concerns the quark-ghost kernel form factor Y_0(q^2) at high energies, the power counting analysis is compatible with 
having a Y_0 (q^2) = 1 at 
large momenta as required by perturbation theory and by the all-orders result summarised in Eq. (<ref>).





The same analysis for the vector component (<ref>) gives, up to logarithmic corrections,

    g_0 (p^2, (p-q)^2 )   Y_0 (q^2)  1/q {⋯}
           +   g_1 (p^2, (p-q)^2 )   Y_1 (q^2)  1/q {⋯}
           +   g_3 (p^2, (p-q)^2 )   Y_3 (q^2)    q  {⋯} ,

where {⋯} stand for finite expressions involving A(p^2), A(q^2), B(p^2) and B(q^2). 
The conditions given in (<ref>) are sufficient to ensure a finite result associate to the UV integration over q for the 
vector component of the Dyson-Schwinger equations. 

The dynamics of QCD generates infrared mass scales for the quark and gluon propagators, see Secs. <ref> and <ref>,
that eliminate possible infinities associated to the low momentum  limit in the integral of the quark gap equation and the analysis of the infrared limit does not 
add any new constraints.

For full QCD, the λ_1 form factor was computed in the soft gluon limit, i.e. vanishing gluon momenta, using 
lattice simulations in <cit.>. The analysis of the lattice data performed in <cit.> shows that the lattice data
is well described by

    λ_1(p^2) = A(p^2) { a + b   D(p^2) }
   
where a and b are constants, that in terms of Y_1 and Y_3 translates into

    2   M(p^2)   Y_1(p^2) - 2   p^2   Y_3 (p^2) ∝ D (p^2)  ,

where M(p^2) = B(p^2) / A(p^2). This result suggests to write

    X_1 ( p^2, (p-q)^2, q^2) = D (  p^2  + (p-q)^2 /2)   Y_1 (q^2)

that in the high q limit gives X_1 ∼ Y_1 (q^2)/q^2 and regularises the ultraviolet behaviour in agreement with
the discussion summarised in (<ref>). 
Similarly, equation (<ref>)  also suggests

    X_3 ( p^2, (p-q)^2, q^2) =  D (  p^2  + (p-q)^2 /2)   Y_3 (q^2)

giving at high q momenta a X_3 ∼X̃_3 (q^2) / q^2 and, in this way, the ultraviolet problems referred in (<ref>) are solved.
Furthermore, for large quark momentum the ansatz (<ref>) and (<ref>) give

    q^3   D (q^2)   F(q^2) /A^2(k^2)   k^2 + B^2(k^2)  p → + ∞   q^3   D (q^2)   F(q^2) /A^2(p^2)   p^2

implying the vanishing of the kernels (<ref>) –  (<ref>) for sufficiently large p. 


In short, our ansatz for the quark-gluon vertex used to solve the Dyson-Schwinger equations reads

    X_0 ( p^2, (p-q)^2, q^2 )     =     X_0( q^2)  ,  
    
       X_1 ( p^2, (p-q)^2, q^2 )     =     D (  p^2 + (p-q)^2/2)    Y_1 (q^2)   , 
    
       X_3 ( p^2, (p-q)^2, q^2 )     =     D (  p^2 + (p-q)^2/2)    Y_3 (q^2)  .

The quark gap equation should be solved taking into account the constraint (<ref>) that demands 

    X_0 (0)    =    X_0 ( q → + ∞ )    =    1  .

The Landau gauge lattice gluon propagator as given by lattice simulations <cit.> retuns a D(q^2) that
is strongly enhanced at low momenta. It follows from Eqs. (<ref>) and (<ref>) that within the ansatz considered here, one expects
X_1 and X_3 to rise significantly for small p^2 + (p-q)^2 = 2   p^2 + q^2 - 2 p · q. On the other hand, at high momenta, the form factors should approach 
its perturbative value. At tree level in perturbation theory the quark-ghost kernel form factors read X_0 = 1 and X_1 = X_3 = 0, suggesting that X_1 and X_3
give  marginal contributions to the full vertex at sufficient high energy. 
At the qualitative level, the guessed behaviour associated with the ansatz (<ref>) – (<ref>) reproduce the computed
quark-ghost kernel form factors computed in <cit.> using the Dyson-Schwinger equations and one-loop dressed perturbation theory
for the quark-ghost kernel.
For the inversion of the Dyson-Schwinger equations, i.e. from the numerical point of view,
given the strong enhancement of the gluon propagator at low momenta, this can mean a poorer resolution of X_1 and X_3 in the ultraviolet regime.




§ PREPARING TO SOLVE THE EUCLIDEAN DYSON-SCHWINGER EQUATIONS 


The computation of a solution of the Dyson-Schwinger equations requires parameterising either the quark-gluon vertex, if one aims to look at the quark propagator,
or the quark propagator functions to extract information on the quark-gluon vertex. In both these cases, a complete description of the
gluon and ghost propagators is assumed explicitly. 

In the current work, we aim to solve the gap equation for 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 fitting the Landau gauge lattice propagators with model functions that are compatible with the results of 1-loop renormalisation group improved 
perturbation theory. 
In this way, it is ensured that the perturbative tails are taken into account properly in the parameterisation of the propagators.
The parameterisations considered here are compared to those of <cit.> in <ref>.
As can be seen on Fig. <ref>, the differences between the two sets of curves are more quantitative than qualitative.





 §.§ 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. The gluon propagator
is well known for the pure Yang-Mills theory and it was calculated in <cit.> for large statistical ensembles and for large physical volumes
∼ (6.6 )^4 and ∼ (8.2 )^4. Furthermore, in <cit.> 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 given there, we will use to solve 
the integral Dyson-Schwinger equations the following fit to the (6.6 )^4 volume result

    D (p^2) = Z   p^2 + M^2_1/p^4+ M^2_2   p^2 + M^4_3[ ω ln(p^2 + m^2_0/Λ^2_QCD) + 1 ]^ γ ,

with the gluon anomalous dimension being
γ = -13/22, Z  = 1.36486 ± 0.00097, M^2_1   =  2.510 ± 0.030 ^2, M^2_2   =  0.471 ± 0.014 ^2, 
M^4_3   =  0.3621 ± 0.0038  ^4, m^2_0   =  0.216 ± 0.026  ^2 using Λ_QCD = 0.425 and where
ω = 33  α_s ( μ ) / 12 π with a strong coupling constant α_s ( μ = 3  ) = 0.3837;  see <cit.> for 
details. This  fit to the lattice data has an associated χ^2/d.o.f. = 3.15. The authors provide fits with better values for the
χ^2/d.o.f. However, given that the level of precision achieved on lattice simulations for the quark propagator is considerably smaller
than for the gluon propagator, one should not distinguish between the various fitting functions provided in <cit.>. Our option
considers the simplest functional form given in that work.

The lattice data for the Landau gauge gluon dressing function p^2 D (p^2), renormalised in the MOM-scheme
at the mass scale μ = 3 GeV and the fit associated to Eq. (<ref>)
can be seen on the top part of Fig. <ref>.

For the ghost propagator we take the data reported in <cit.> for the 80^4 lattice simulation and fit the lattice data to the functional
form

    D_gh (p^2)  =  F(p^2)/p^2
           =  Z/p^2  p^4 + M^2_2   p^2 + M^4_1/p^4+ M^2_4   p^2 + M^4_3 [ ω ln(p^2 + m^4_1/p^2 + m^2_0/Λ^2_QCD) + 1 ]^ γ_gh ,

getting Z  =  1.0429 ± 0.0054, M^4_1 = 18.2 ± 5.7 ^4, M^2_2 = 33.4 ± 6.4 ^2,
M^4_3 = 6.0 ± 2.7 ^4, M^2_4 = 29.5 ± 5.7 GeV^2,
m^4_1 = 0.237 ± 0.049, m^2_0 = 0.09 ± 0.42 ^2
with a χ^2/d.o.f. = 0.27. In the above expression the ghost anomalous dimension reads γ_gh = - 9/44 with
ω and Λ_QCD taking the same values as in the gluon fitting function (<ref>).
The lattice data, renormalised in the MOM-scheme at the mass scale μ = 3 GeV,  and the fitting curve (<ref>) can be
seen on the bottom of Fig. <ref>.



 

 §.§ Lattice Quark Propagator 


For the quark propagator we consider the result of a N_f = 2 full QCD simulation in the Landau gauge <cit.> 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  <cit.> and  <cit.>
for 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 <cit.>, 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. <ref>. The UV corrected lattice data
is then fitted to the rational function

    Z(p^2) = Z_0   p^4 + M^2_2   p^2 + M^4_1 /p^4 + M^2_4   p^2 + M^4_3
 
giving Z_0 = 1.11824 ± 0.00036, M^4_1 = 1.41 ± 0.18 ^4, 
M^2_2 = 6.28 ±1.00 ^2, M^4_3 = 2.11 ± 0.28 ^4,
M^2_4 = 6.20 ± 0.98 GeV^2 for a χ^2/d.o.f. = 0.74. The solid red line on Fig. <ref> (top) 
refers to the fit just described.

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 <cit.>. 
The lattice data published in <cit.> and reported on Fig. <ref> (bottom) was obtained using the 
so called hybrid corrections to  reduce the lattice effects <cit.> . 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; see Appendix on <cit.>. 
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. <ref> (bottom) 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. <ref> the data in the region with an orange background was
not taken into account in the global fit of the running quark mass. The remaining lattice data was fitted to

    M(p^2) =  m_q ( p^2) /[ A + log( p^2 + λ  m^2_q(p^2) ) ]^γ_m

where γ_m = 12/29 is the quark anomalous dimension for N_f = 2 and

    m_q(p^2) = M_q   p^2 + m^2_1/ p^4 + m^2_2   p^2 + m^4_3  + m_0  .

The fitted parameters are
 M_q = 349 ± 10 ^2, m^2_1  = 1.09 ± 0.43 ^2, m^2_2  = 0.92 ± 0.28 ^2, 
m^4_3 = 0.42 ± 0.15 ^4,
m_0 = 10.34 ± 0.63 MeV and A = -2.98 ± 0.25 for a χ^2/d.o.f. = 1.97 after setting λ = 1 ^2/^2.
 The fit function and the full lattice running
quark mass data can be seen on Fig. <ref> (bottom). Note that in (<ref>) 
and in (<ref>) p is given in GeV and m_q(p^2) and M(p^2) are given in MeV.





§ SOLVING THE DYSON-SCHWINGER EQUATIONS 




Let us now discuss the solutions of the Euclidean space Dy­son-Schwinger equations (<ref>) and (<ref>) 
for the quark-ghost kernel, i.e. for the quark-gluon vertex. The momentum integration will be performed as described in <ref>, 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 <cit.>.
After angular momentum integration, one is left with the kernels 

    𝒩^(0,1,3)_A,B (p,q) =   q^3  D(q^2)   F(q^2)/[ A(k^2) ]^2    k^2 + [ B(k^2) ]^2 𝒦^(0,1,3)_A,B (p,q)

that can be seen on Figs. <ref>, <ref> and <ref>, 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. <ref>, <ref> and <ref> and the main difference being that the associated numerical values are considerably smaller.






The major contributions of the 𝒩^(0)_B and 𝒩^(0)_A kernels occur
in a well defined momentum region where p ≲ 1.5 GeV and q ≲ 1.5 GeV. Further, for p, q ≳ 2 GeV the kernels become marginal. 
The results for 𝒩^(0)_B and 𝒩^(0)_A  reproduce the corresponding behaviour observed in <cit.>.
It follows that the integration over momentum in Eqs. (<ref>) and (<ref>) associated to the kernels 
𝒩^(0)_B and 𝒩^(0)_A kernels, that  are coupled to X_0(q^2), is finite.

The function 𝒩^(1)_A (p, q) displays a similar pattern and, again, the integration over the gluon momentum associated with 
𝒩^(1)_A is expected to be well behaved. 
On the other hand the remaining kernels, i.e.  𝒩^(1)_B (p, q), 𝒩^(3)_B  (p, q) and 𝒩^(3)_A  (p, q),
are all increasing functions of q. The requirement of a finite integration over q demands that X_1 and X_3 should approach zero fast enough
to compensate the increase with q of these kernel functions; see the discussion of the kernels ultraviolet limit in Sec. <ref>. 
The ansatz (<ref>) – (<ref>) adds a multiplicative gluon propagator term that is just enough to regularize the 
high momentum associated to 𝒩^(1)_B (p, q), 𝒩^(3)_B  (p, q) and 𝒩^(3)_A  (p, q).
Indeed if one takes into account the multiplicative gluon propagator contribution to the kernels, those who are divergent become well behaved.
This can be seen on Figs. <ref> and <ref> where the kernels, now including the multiplicative gluon propagator term,
are reported. The new versions of 𝒩^(1)_B (p, q), 𝒩^(3)_B  (p, q) and 𝒩^(3)_A  (p, q)
mimic the pattern observed for 𝒩^(0)_B, 𝒩^(0)_A and 𝒩^(1)_A (p, q) and, once more, their main contribution to the integral
equations happens for p ≲ 2 GeV and q ≲ 2 GeV. 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. All quantities are renormalised in the MOM scheme,
using the same renormalisation scale as in <cit.>, i.e μ = 4.3 GeV, so that one can compare easily the results of the two works.
The renormalized quantites satisfy the identities

    Z ( μ^2 )  =  1/A( μ^2 ) = 1  ,    D(μ^2) =  1/μ^2 ,    F( μ^2 ) = 1  .

The bare quark mass quoted in the lattice simulation for the ensemble used here reads m^bm = 8 MeV <cit.>.
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 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. <ref>, <ref>, <ref> and <ref> were computed using the same
value for α_s( μ ) = 0.295 as in <cit.>. In Sec. <ref> we allow α_s( μ ) to deviate from this value
and provide a “best value”. From Sec. <ref> onwards, the results reported use the optimal value for the strong coupling constant.








 §.§ One-Loop Dressed Perturbation Theory for X_0(q^2) 





The four longitudinal quark-gluon form factors were pa­ra­me­tri­sed in terms of the three quark-ghost kernel form factors X_0, Y_1 and Y_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 one sets Y_1 = Y_3 = 0 and, then, solve the gap equation to estimate Y_1 and Y_3. The way 
the solutions of the Dyson-Schwinger equations for Y_1 and Y_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. <ref> that, in the simplified version of
kernel, translates into the following integral equation

    X_0(p^2)  =  1 
           + C_F   g^2/8∫_k^Λ [ p^2 - (k · p)^2/p^2]   D(k^2)   F ( (k + p/2)^2 )   F(k^2)/(k + p/2)^2×
          ×  A( (k+p)^2 )  [  A( (k+p)^2 ) + A(p^2) ]   H_1( (k+p)^2 ) /A^2( (k+p)^2 )   (k+p)^2  + B^2( (k+p)^2 ) ×
          × 
        X_0 (k^2)

with  H_1 (q^2) representing 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 <cit.> where

    H_1(q^2) =  c ( 1 +  a^2 q^2/q^4 + b^4) + ( 1 - c ) w^4/w^4 + q^4 ,
  
for c = 1.26, a = 0.80 GeV, b = 1.3 GeV and w = 0.65 GeV. 
For a recent analysis of the quark-ghost vertex see <cit.>.
Equation (<ref>) was solved after introducing a cutoff Λ = 100 GeV, after performing
the angular integration using 1000 Gauss-Legendre points and considering 2000 Gauss-Legendre points for the integration over k. The introduction of a
Gauss-Legendre quadrature reduces the integral equation (<ref>) 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. <ref> and are, essentially, those reported in <cit.>.
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 Eq. (<ref>) but no convergence was observed, therefore, the solutions reported  on
Fig. <ref> are those computed from a single iteration.



The estimation of X_0 allows to solve the gap equation for Y_1 and Y_3. In order to solve the Dyson-Schwinger equations, after
the angular integration, the scalar and vector components of the equations are rewritten in the form of the larger linear system

    ( [ B(p) - Z_2 m^b.m. -  𝒩^(0)_B X_0 (p);                                     ;       A(p) - Z_2 -  ⁄𝒩^(0)_A X_0 (p) ]) = 
           = 
    ( [ 𝒩^(1)_B 𝒩^(3)_B;                ; 𝒩^(1)_A 𝒩^(3)_A ])  ( [ Y_1;    ; Y_3 ])    ;

from now on we will adopt the short name version B = 𝒩 X to refer to this linear system of equations. Note that 
Y_1 has mass dimensions, while X_0 and Y_3 are dimensionless. Note that the kernels in 𝒩
also have different dimensions and it is only after multiplication that we recover the proper dimensionful equation.

A direct solution of B = 𝒩 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 𝒩, that translates the ill defined problem in  hands. 
The linear system can be solved using the Tikhonov regularisation <cit.> that  replaces the original linear system by a minimisation 
of the functional || B - 𝒩 X ||^2 + ϵ ||X||^2, where ϵ is a small parameter to be determined in the inversion. 
This functional favours solutions that solve approximately the linear system but whose norm is small. 
For real symmetric matrices, Tikhonov regularisation replaces the original
system by  its normal form 
𝒩^T B = ( 𝒩^T 𝒩  + ϵ )X. 
Although in our case 𝒩 is not a symmetric matrix, we will solve the system as given in its normal form.



The determination of the optimal ϵ is done by solving 𝒩^T B = ( 𝒩^T 𝒩  + ϵ )X for various values of
ϵ and look at how || B - 𝒩 X ||^2 and ||X||^2 behave as a function of the regularisation parameter ϵ. 
The outcome of the inversions for different ϵ
can be seen on Fig. <ref>. 
For smaller values of ϵ, i.e. when one is closer to the original ill defined problem, the corresponding
solution of the linear system results on Y_1 and Y_3 with larger norms.
The larger values of the regularisation parameter ϵ are associated to solutions of the modified linear system with smaller Y_1 and Y_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 equations,
is among the smallest values just before the norms of Y_1 and Y_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. <ref> 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 Y_1(p^2) and Y_3(p^2)
from solving the modified gap equations, provided we let the norm of Y_1 and Y_3 to be large enough. Of course, for large norms Y_1 and Y_3
are free to vary over a large range of values and the solutions with smaller norms are preferred.

From Fig. <ref>  three typical solutions 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




    
                   ||Y_1||^2      ||Y_3||^2      ||Δ||^2      ||Δ||^2      ϵ   
 
   I (TL)         1.991              5.945              0.00125                               0.03811                                0.0095         

   II (TL)        0.999              0.594               0.01739                                0.3151                                0.175           

   III (TL)        0.749           0.324               0.03796                                0.4113                                 0.291           
   
   



for m^b.m. = 6.852 MeV, Z_2 = 1.0016, where Δ = B - Z_2 m^b.m. - K_B X, ||Y_2|| is given in GeV,
and Δ = A - Z_2 - K_A X with ||Y_3|| being dimensionless. For all these solutions ||X_0 - 1||^2 = 0.24214.
On the other hand, the characteristics of the solutions computed with the enhanced (Enh) ghost-gluon vertex given by Eq. (<ref>) are




    
                        ||Y_1||^2      ||Y_3||^2      ||Δ||^2      ||Δ||^2      ϵ   
 
   I (Enh)           1.994             6.464               0.001257                             0.03526                               0.0085           

   II (Enh)          0.999             0.633               0.01644                               0.3178                                0.17          

   III (Enh)         0.745            0.350               0.03536                               0.4146                                0.28           
   
   



have the same m^b.m. and Z_2 as the previous ones and ||X_0 - 1||^2 = 0.45283. In both cases, the norms of X_0, Y_1 and Y_3 are the norms of the 
corresponding part of the vector that appear in the linear system .

The quality of the solutions can be appreciated on
Fig. <ref> 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 Y_1 and Y_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. <ref>. On the figures we have defined

    Δ    =     B - (Z_2 m^b.m. + K^(0)_B X_0 + K^(1)_B Y_1 + K^(3)_B Y_3 )  ,  
    Δ    =     A - (Z_2 + K^(0)_A X_0 + K^(1)_A Y_1 + K^(3)_A Y_3 )  .

||Δ||^2  and ||Δ||^2  should be understood as the sum of the squares of the components of the linear systems
(<ref>) and (<ref>), respectively, over the Gauss-Legendre points.
As Fig. <ref> 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 ||Δ||^2 relative to ||Δ||^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 Y_1 and Y_3 that don't seem reasonable and, for example, result in
a Y_1 at the cutoff that is far away from zero. Further, for the rescaled systems the Δ and Δ 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 Y_1 and Y_3 computed for the various ϵ associated to the solutions I (TL) – III (TL) and I (Enh) – III (Enh)
can be seen on Fig <ref>. For Y_1 the outcome of resolving the integral equations using either the tree level or the enhanced
ghost-gluon vertex result on essentially the same function. Further, the various solutions provide essentially the same Y_1(p^2), with the exception of
III (TL) and III (Enh) that return a suppressed form factor relative to the other solutions. For Y_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. <ref>, 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%.

The quark-ghost kernel
form factors X_0, X_1 and X_3 were also computed in <cit.>, see their Fig. 4, combining the quark gap equation with a one-loop
dressed perturbation theory for the quark-ghost kernel. The comparison between the results of the two calculations is not straightforward. 
Indeed, if in our calculation X_0 is assumed to be a function only of the gluon momentum, in <cit.> the authors take into account 
its full momentum dependence and evaluate how it changes with the quark momentum, the gluon momentum and the angle between these two momenta. 
There calculation results on a X_0 that is always close to its tree level value X_0 = 1 and whose values are in the range [1 , 1.12 ].
A qualitative comparison with the here reported form factor, shows that, from the point of view of the dynamical range of values, our X_0 computed using the 
enhanced ghost-gluon vertex is closer to that reported in <cit.>. In both cases X_0 is always close to its tree level value and differs from unit by, at most,
10%.

The comparison between the remaining form factors is slightly more involved. Indeed, the X_1, X_2 and X_3 referred in <cit.>
compared with the expressions given in Eqs. (<ref>) and (<ref>)  and not directly with the Y_1 and Y_3 reported in 
Fig. <ref>. Note that there are signs differences on the definition of the various quark-ghost kernel form factors between <cit.>
and the current work.
Due to the presence of the gluon propagator in (<ref>) and (<ref>) one expects some angular dependence of the quark-ghost kernel form factors.
Further, due to the parameterisation used for the gluon propagator, the form factors are expected to be enhanced when, simultaneous, the quark and gluon momenta 
become smaller, i.e. in the infrared limit.
This is precisely what is observed for the X_1 and X_2 reported in <cit.>.




Our estimations for X_1 and X_2 for the solutions referred previously and for the particular kinematics p = 0, sometimes called the soft quark limit, can be seen 
on Fig. <ref>. The first comment about these form factors is that they seem to be independent of the ghost-gluon vertex and, indeed,
the form factors associated to the solution using a tree level ghost-gluon vertex, named (TL) in the figure, are essentially indistinguishable from those computed
using the enhanced ghost-gluon vertex, named (Enh) in the figure.
The quark-ghost kernels form factors X_1 and X_3 differ from there perturbative values at low momenta, i.e. for q ≲ 1 GeV for X_1 and
for q ≲ 2 GeV for X_3. According to <cit.>, these two form factors increase (in absolute value) for momenta
below ∼ 1 GeV, a result that is in good qualitative agreement with our estimations. 
From their Fig. 4, it is not clear if X_3  0 extends over a wider range of momenta, when compared to X_1. Our calculation returns a X_3 that differ from zero on
larger range of momenta, when compared with X_1. 
Furthermore, the X_1 and X_3 computed in <cit.> are monotonic increase functions (absolute values) when the zero momentum
limit is approached. The form factors reported on Fig. <ref> show a pattern of maxima, with X_1 having a single maximum for q ∼ 350 MeV for X_1 and
X_3 showing several maxima (in absolute value) at momenta q ∼ 250 MeV, ∼ 600 MeV and ∼ 1 GeV. Note that the zero crossing for X_3 occur
for momenta that are of the same order of magnitude of Λ_QCD, the gluon mass (or twice Λ_QCD) and the usual considered as a non-perturbative
mass scale (1 GeV). It is not obvious why the zero crossing of X_3 occur for such mass scales and why the crossing is not
seem for X_1. If the form factors reported on Fig. 4 of <cit.> never cross the zero value, that is not the case of the form factors represented on Fig  <ref>. X_1 shows various zeros that, curiously, seem to disappear for the solution with the smaller norm.
On the other hand, X_3 cross  zero  for q ≈ 363 MeV and 835 MeV for solutions II and III. This is a major difference between the two sets of
solutions under discussion. Another important difference being the dynamical range of values. Our X_1 is within the range of values [ 0 , 3.5] GeV^-1, while the same
form factor computed in <cit.> is within [ 0 , 0.2] GeV^-1 which represents a factor of ∼ 20 smaller than our estimation.
However, the two calculations report a |X_3| within the range [ 0 , 0.5] GeV^-2. Our result overestimates X_1 relative to <cit.> but
returns a |X_3| within the same dynamical range of values. Another major difference being that the maxima of the form factors does not occur at vanishing momenta
as in <cit.> but at finite and small momenta, ∼Λ_QCD for X_1 and ∼ 2  Λ_QCD for the absolute maxima of |X_3|.


The estimations of X_0, X_1 and X_3 suggest that the quark-gluon vertex is dominated, at the infrared by those terms that are associated with
X_1. If this is the case, then, given the definitions (<ref>) – (<ref>) one expects that the dominant contributions to the quark-gluon vertex
to be associated with the form factors

    λ_1 (-p,   p-q,   q )  ≅  F(q^2)   [ B(p^2)  + B( k^2 ) ] ×
          ×  D( p^2 + (p-q)^2/2)   Y_1(q^2) ,  
    
          λ_3 (-p,   p-q,   q )  ≅
          ≅  F(q^2)   A( k^2 )  ( q^2 -   (p q) ) - A(p^2)   ( p q ) /q^2 - 2 (p · q)×
          × D( p^2 + (p-q)^2/2)   Y_1(q^2)  , 
    
          λ_4 (-p,   p-q,   q ) ≅
          ≅F(q^2)/2  Δ A    D( p^2 + (p-q)^2/2)    Y_1(q^2)

and since Z ∼ 1 and B ∝ M(p^2), one expects the dominant form factor for the quark-gluon vertex to be associated with the tree level operator, i.e. with λ_1.  






 §.§ 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 regularisation of the
corresponding linear system of equations follow the steps described in the previous section.
First the angular integration is performed using 800 Gauss-Legendre points. 
Then, for the momentum integration the cutoff Λ = 20 GeV is introduced and  we consider 200 Gauss-Legendre points to perform 
the integration over the loop momentum. Further, to determine the solutions of the Dyson-Schwinger equations, now already in the form of a 
linear system of equations, the scalar component and the vector component of the gap equation are grouped into a large linear system as follows

    ( [ B(p) - Z_2 m^b.m.;                  ;        A(p) - Z_2 ]) =
    ( [ 𝒩^(0)_B 𝒩^(1)_B 𝒩^(3)_B;                ; 𝒩^(0)_A 𝒩^(1)_A 𝒩^(3)_A ]) 
    ( [ X_0; Y_1; Y_3 ])

that again we refer, as a short name, by B = 𝒩 X. The upper  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 Y_1 and Y_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 (<ref>) returns
X_0(q^2) for q ∈ [ 0 ,  Λ] and Y_1(q^2) and Y_3(q^2) for q ∈ [ 0 ,  Λ / 2]. 
For Y_1 and Y_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 for X̃_0(q^2), 
rebuilding X_0 (q^2) at the end.
The resulting linear system is then regularised using the Tikhonov regularisation  and the corresponding 
𝒩^T B = ( 𝒩^T 𝒩  + ϵ )X linear system is solved for various ϵ.
The choice of the optimal regularisation parameter ϵ follows the criteria discussed in Sec. <ref>. 
We have checked that by interchanging the roles of X_0, Y_1 and Y_3 when building the large linear system the solutions are unchanged;
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. <ref> the residuum squared of the scalar (top) and vector (bottom) components of the gap equation are shown against the norm
of various form factors. Smaller values of ϵ are associated to solutions with larger norms and appear at the right side of the plots, while larger values
of ϵ are associated to solutions with smaller norms that show up on the left side of the plots. As shown on the figure, the residuum for both equations
has a stronger dependence on ϵ and can take quite small values. For the solutions featured in the plot, the smallest residuum squared reaches values of
the order of 10^-5 for the scalar equation and 10^-4 for its vector component. Similar values for the minimum residuum were also observed for the 
solutions computed using the perturbative estimation of X_0.


On Fig. <ref> we identify four solutions associated to an ϵ around its optimal value and whose characteristics are




    
                    ||X_0 - 1||^2       ||Y_1||^2       ||Y_3||^2      ||Δ||^2      ||Δ||^2      ϵ   
 
   I                0.5003                   2.0039              4.2966             0.001205                               0.05568                               0.012   

  II                0.09298                 1.4994              1.2538              0.004890                               0.1730                                 0.071   

  III               0.03408                 1.0002              0.4634             0.02034                                 0.3154                                 0.191   

  IV               0.01811                 0.7001              0.2150             0.05486                                 0.4291                                 0.365   

   
   



for m^b.m. = 6.852 MeV, Z_2 = 1.0016, where Δ and ||Y_1|| are  given in GeV, while Δ, ||X_0|| and ||Y_3||
are dimensionless.




The relative errors for the solutions I – IV of the regularised linear system are show on Fig. <ref>.
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
non-perturbative 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 in Sec. <ref>. The observed growth of the relative error for p ≳ 10 GeV is probably related also to the
missing components of Y_1 and Y_3 which are set to zero for this range of momenta.


The form factors X_0, Y_1 and Y_3 associated to the solutions I – IV can be seen on Figs. <ref> – <ref>.
On Fig. <ref> besides solutions I –IV we also show the perturbative X_0(q^2) computed using one-loop dressed perturbation
theory with the tree level ghost-gluon 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 regularised 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 non-perturbative 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  resolve the gap equation with the smaller relative errors that is below 8%
– see Fig. <ref>.
 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 take 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 Y_1(p^2) can be view­ed on Fig. <ref>. All the solutions I – IV reproduce the same pattern for
this form factor, with a positive maxima around p ≃ 400 MeV and with Y_1 becoming small for p ≳ 1.5 GeV. In particular, for the solutions
III and IV, Y_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 Y_1 but this function times the gluon propagator – see Eq. (<ref>). The same applies to
Y_3 as can be seen on  Eq. (<ref>). Once more, as the norm of Y_1 decreases, the form factors seems to prefer to take only positive values.



The form factor Y_3(p^2) is reported on Fig. <ref>. It turns out that this function is positive for p ≲ 400 MeV and for p ≳ 1 GeV,
takes negative values in between, with a maximum at p ≃ 1.5 GeV, and then slowly approaches its tree level value from above. Given the way
the solutions are computed, on Fig. <ref> Y_3(p^2) shows a jump at p ≃ 10 GeV that corresponds to p = Λ/2. A similar
behaviour can be seen on Fig. <ref> for Y_1(p^2). However, given that for p ≃ 10 GeV  one has a Y_1(p^2) ≃ 0, this sudden jump is not
so easily observed.



Finally, on Fig. <ref> we provide the various solutions for X_0(p^2), Y_1(p^2) and Y_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 Y_1 and Y_3 being obtained only in the range p ∈ [ 0   ,  Λ/2].
For the so-called X_1   X_0   X_3 the extended vector included Y_1 over the full set
of Gauss-Legendre points with X_0 and Y_3 being obtained only in the range p ∈ [ 0   ,  Λ/2].
For the so-called X_3   X_0   X_1 the extended vector included Y_3 over the full set
of Gauss-Legendre points with X_0 and Y_1 being obtained only in the range p ∈ [ 0   ,  Λ/2]. 
We call the readers attention to the stability of the solution of the various linear
systems. Furthermore, the comparison of Fig. <ref> from Sec.<ref> and Fig. <ref>
show quite similar Y_1 and Y_3 suggesting, once again, that X_0 almost does not deviates from its tree level value.





 §.§ Solving the DSE for X_0 = 1 








The non-perturbative solutions of the Dyson-Schwinger equations discussed on the previous paragraph suggest that X_0(p^2) ≃ 1. Therefore,
herein we investigate the results by solving the DSE with X_0(p^2) = 1. The residua of the scalar and vector equations against the norm of the two remaining form factor Y_1 and Y_3 can be seen on Fig. <ref>. The characteristics of the solutions
highlighted in the figure and associated to an ϵ close to its optimal value are




    
                   ||Y_1||^2      ||Y_3||^2      ||Δ||^2      ||Δ||^2      ϵ   
 
   I               2.4870              10.0286            0.0004739                          0.01314                                0.0025   

  II               2.0003                4.6676         0.001351                             0.04691                               0.013   

  III              1.5006                1.3225         0.005306                             0.1695                                  0.0745   

  IV              1.0009                0.4801         0.02232                               0.3157                                  0.1975   

   
   



for m^b.m. = 6.852 MeV, Z_2 = 1.0016, where Δ and ||Y_1|| are  given in GeV, while Δ, ||X_0|| and ||Y_3||
are dimensionless. The relative error on the Dyson-Schwinger equations for these solutions can be seen on Fig. <ref> 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 Y_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 Y_1(p^2) and Y_3(p^2) associated to the solutions I – IV are reported on Fig. <ref> and reproduce the same
patterns as  the solutions computed in the previous sections.





 §.§ Full Form Factors and Comparison of Solutions 


In Secs. <ref>, <ref> and <ref> we have solved the Dyson-Schwinger equations assuming that the quark-ghost
kernel form factors are given by Eqs. (<ref>) – (<ref>). So, besides, X_0, the full form factors appearing in H and H,
see Eqs. (<ref>), 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. <ref>. This figure repeats 
partially Fig. <ref> 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, with the perturbative solutions having the largest deviation from unit.



The form factor Y_1(p^2) can be seen on Fig. <ref> for all the solutions. Note that all solutions reproduce essentially the same
function of the gluon momentum, with Y_1(p^2) being small for p≳ 1.5 GeV and showing a sharp peak at p ≃ 400 MeV. Y_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 Y_3(p^2) can be seen on Fig. <ref> for all the solutions. Surprisingly, Y_3(p^2) seems to have a relative large tail that
appears in 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 non-perturbative solution obtained from inverting the gap equation. Y_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 negative values in p ∈ [ 0.35   ,   0.9] GeV.

In summary, Figs. <ref> – <ref> resume the computations of the quark-ghost kernel form factors performed so far.





 §.§ Tunning α_s 


The results for the relative errors on the scalar and vector components of the Dyson-Schwinger equation 
seen on Figs. <ref>, <ref>, <ref> 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 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, (<ref>), is set to identity[It can also be viewed as been
included in the definition of the various form factors.]. 
Although the original integral equation is linear on the form factors X_0, Y_1 and Y_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 regularised linear system of equations that replace the original equations
in the way it is described on Sec. <ref> for α_s(μ) = 0.20, 0.22 and 0.25. The properties of the inversions of the regularised
system for the various values of the strong coupling constant can be found in Figs. <ref>,
<ref> and <ref> and should be compared with Fig. <ref> of the Sec. <ref>.







As the figures shows, lowering the value of α_s( μ ) solves the problem of the increase of the relative error observed in Sec. <ref>. 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. <ref>. The relative error associated to
the remaining solutions given on Figs. <ref>, <ref>, <ref> and 
<ref> are larger and, therefore, we take α_s(μ) = 0.22 as the optimal value for the strong coupling constant
within our approach. The corresponding quark-ghost kernel can be seen on Figs. <ref>, <ref> and <ref>, 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), Y_1(q^2), Y_3(q^2) that, together with
the quark, gluon and ghost propagators,  define the full form factors as given in Eqs. (<ref>),  (<ref>) and (<ref>).
Once the full quark-ghost kernel form factors are known, then the longitudinal quark-gluon form factors can be computed using
Eqs. (<ref>) – (<ref>), 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         
            g_1(p^2_1, p^2_2) = g_2(p^2_1, p^2_2) = D (  p^2_1 + p^2_2/2)  .

For completeness, we write the full expressions for the longitudinal form factors in Euclidean space

    λ_1 (-p,   k = p-q,   q  ) =  F(q^2)/2{
          [ A(p^2)  + A( k^2 ) ]    X_0(q^2) 
           
                  +   2  [ B(p^2)  + B( k^2 ) ]   D( p^2 + k^2/2)   Y_1(q^2)  
           +  [ A(p^2) ( (pq) - 2 p^2 )  
           + A(k^2) ( 3 (pq) - 2 p^2 - q^2 )   ]  ×
          ×  D( p^2 + k^2/2)   Y_3(q^2)  }  ,


    λ_2 (-p,   k = p-q,   q )  =  F(q^2)/2( 2 (p · q) - q^2 ){
          [ A(k^2) - A(p^2) ]    X_0(q^2)  
                  + [  A(k^2) ( q^2 -  (pq) ) - A(p^2)   (pq)   ]  ×
          × D( p^2 + k^2/2)   Y_3(q^2) }  ,


    λ_3 (-p,   k = p-q,   q ) =  F(q^2)/2 (p · q) - q^2 {
          [ B(k^2) - B(p^2) ]      X_0(q^2) 
           + [  A(p^2)   ( p q )  - A( k^2 )  ( q^2 -   (p q) )   ] ×
          × D( p^2 + k^2/2)    Y_1(q^2) }  ,


    λ_4 (-p,   k = p-q,   q ) =  F(q^2)/2   D( p^2 + k^2/2)  ×
          ×{[ A(k^2) - A(p^2) ]      Y_1(q^2)  
          + [ B(k^2) - B(p^2) ]    Y_3(q^2) }  .

Note that by taking into account  structures of the quark-ghost kernel other than X_0 the quark-gluon vertex deviates considerably 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 D( (p^2 + k^2)/2 ).

For the calculation of λ_1–λ_4 we will use Sol. II computed using α_s( μ ) = 0.22; see Sec. <ref> 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. <ref>. On Fig. <ref> 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^(1)_μ – L^(4)_μ, see Eqs. (<ref>) for their definition in Minkowski 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 a minor role, i.e. the tensor structures of the longitudinal
part of the vertex seem to be sub-leading; 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 <cit.>, where all the strength appears 
associated to λ_1. 
For example, for the kinematical configuration defined by p^2 = (p-q)^2 at vanishing p they 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 renormalisation 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. <ref> and <ref> 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 <cit.> is difficult to perform but in these 
works λ_1 clearly dominates. On <cit.> λ_2 reaches at most 16% of the maximum value of λ_1, while
λ_3 seems to have the possibility of taking large values. On <cit.>, λ_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. <ref>.

As seen on Figs. <ref> and <ref> 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 (<ref>) to
rewrite the form factor. Indeed, as can be seen on Eqs. (<ref>) – (<ref>), all the form factors have, as a global factor, the ghost
dressing function F(q^2). The ghost dressing function is enhanced, roughly by a factor of three, 
in the infrared, see Fig. <ref>,  implying the increase of the λ_i as q = 0 is approached.

The infrared enhanced of the form factors with the quark incoming momentum is more subtle. It is linked to our ansatz 
that relies on the analysis of the soft gluon limit of the Landau gauge lattice data for λ_1 performed in <cit.>.
Indeed, this work identified a dependence of λ_1 on the gluon propagator that was incorporated in the ansatz,
making the quark-ghost kernel form factors X_1 and X_3  proportional to D ( (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 together with the increase of the ghost dressing function.











The dependence of the quark-gluon form factors in the angle between p and q 
can be seen on Figs. <ref> – <ref>. 
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[For example, for p = 0.5 GeV λ_1 and λ_3 
there is an enhancement of a factor of ∼ 2.7 and ∼ 2.0, respectively, between the maxima values for θ = 0 relative to θ = π.
For λ_2 and λ_4 this enhancement is ∼ 3.5 and ∼ 3.2. The corresponding factors for an incoming quark momentum
p = 1 GeV are ∼ 1.5 for λ_1, ∼ 4.4 for λ_2, ∼ 2.3 for λ_3 and a suppression by a factor of ∼ 0.6
for λ_4. Also, for λ_1 and λ_3 the maxima at θ = π/2 is about half of the maxima at θ = 0.].
It follows that the quark-gluon vertex favours small quark and gluon momenta and parallel four-vectors p and q.

From the point of view of the momentum dependence, our solution for the quark-gluon vertex is closer to that of the Maris-Tandy model <cit.>
than those computed in <cit.>, in the sense that we observe a rather strong
enhancement at low momenta. Indeed, compared to these last references, the herein computed form factors are significantly larger.
Recall that the Maris-Tandy model considers a single form factor, that
would be (effective) 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 our quark-gluon form factors are enhanced occurs essentially within the same range of momenta
as  the corresponding effective 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 their maximum values – see Figs. <ref>, <ref>
and <ref>.




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. <ref> (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,
as observed also when studying the solutions associated with the perturbative X_0 as discussed at the end of Sec. <ref>,
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 has recently being computed using full QCD lattice simulations <cit.>.
The data was investigated in <cit.> 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. <ref> – see the full curve in black.
This curve was (arbitrarely) normalised to reproduce the lattice data at 1 GeV of the β = 5.29 and
M_π = 295 MeV simulation[The normalisation essentially removes the F(0) that can appear in the expression for λ_1; see Eq. (<ref>).]. 
Clearly, our ansatz underestimates λ_1 in the infrared region. As discussed in <cit.>, in the soft gluon
limit

    λ_1(p^2) = F(0)/Z(p^2){ 1 + 2   M(p^2)   X_1(p^2) - 2   p^2   X_3(p^2) } ,

where p is the quark incoming momenta, and in our notation

    X_1(p^2) → D ( p^2 )   Y_1(0)
          
      X_3(p^2) → D( p^2 )   Y_3(0)  .

Our solutions has Y_1(0) ≈ 0 GeV and Y_3 (0) ≈ 0 and, therefore, it underestimates λ_1 (p^2)
in the infrared region.
Note that herein Y_1 and Y_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 that, possibly, prevent the inversion to resolve correctly Y_1(q^2) and Y_3(q^2) in the deep infrared region. 
If in the calculation of the soft gluon limit one assumes that Y_1(0) deviates from zero by a small quantity, the agreement with the lattice data
is considerably improved both in the infrared and in the ultraviolet. This is represented by the two full curves in colour of Fig. <ref> where Y_1(0) is
set to a small value. The colour curves suggest a Y_1(0) ∼ 0.05 – 0.07 GeV. 
Further, the agreement in the ultraviolet region can also be improved if Y_3(0) assumes small and positive values; recall that Y_3(q^2) approaches 
zero from the above when q^2 → 0 as can be seen on Fig. <ref>.






§ SUMMARY AND CONCLUSIONS 


In this work we investigated 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 derived in <cit.>. 
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 <cit.>.
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 <cit.>. 
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 <cit.>
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 needs to be regularised 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 regularisation that is equivalent to minimize || B - 𝒩 X ||^2 + ϵ ||X||^2. 
The solutions are found numerically after writing the regularised 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 - 𝒩 X ||^2, and the norm of the corresponding quark-ghost form factors, i.e. ||X||^2, for each solution of the regularised 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. (<ref>),
(<ref>) and (<ref>) – (<ref>) 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 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
 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 effective vertex introduced in <cit.>, and those based on first principles ab initio
continuum methods, see e.g. <cit.> 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. <cit.> and references therein.
  




§ ACKNWLODGEMENTS


This work was partly supported by the
Fundação de Amparo à Pesquisa do Estado de São Paulo [FAPESP Grant No. 17/05660-0], Conselho Nacional de Desenvolvimento Científico e 
Tecnológico [CNPq Grant 308486/2015-3] and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) of Brazil. 
This work is a part of the project INCT-FNA Proc. No. 464898/2014-5. 
O. Oliveira acknowledge support from CAPES process 88887.156738/2017-00 and FAPESP Grant Number 2017/01142-4. 
JPCBM acknowledges CNPq, contracts 401322/2014-9 and 308025/2015-6. This work was part of the project FAPESP Temático, No. 2017/05660-0.
 




§ 4D SPHERICAL COORDINATES AND INTEGRATION OVER MOMENTUM 


In 4D the spherical coordinates are related to the cartesian coordinates as follows

    x_1     =     r  cosϕ_1 
    
     x_2     =     r  sinϕ_1  cosϕ_2 
    
     x_3     =     r  sinϕ_1  sinϕ_2  cosϕ_3 
    
     x_4     =     r  sinϕ_1  sinϕ_2  sinϕ_3

where

    r ∈ [0 ,   + ∞ [   ,     ϕ_1,  ϕ_2 ∈ [0 ,  π]         ϕ_3 ∈ [ 0 ,   2 π [   .
 
The 4D volume element reads

    dV = r^3 dr  ( sin^2 ϕ_1   d ϕ_1 )  ( sinϕ_2   d ϕ_2 )  ( d ϕ_3 )  .
 
Setting the outgoing quark momenta p = (p,   0,   0 ,   0 ) it follows that

    p · q = p   q  cosϕ_1

and the angular integration in the Dyson-Schwinger equations can be written as

    ∫^π_0 sin^2 ϕ_1   d ϕ_1  ∫^π_0 sinϕ_2   d ϕ_2  ∫^2 π_0 d ϕ_3     =     4 π ∫^1_-1sinϕ_1   d ( cosϕ_1 )
      
        =     4 π ∫^1_-1√( 1 - x^2)  dx

where, in the last identity, we set x = cosϕ_1. It follows that momentum integration in the Dyson-Schwinger equations reads

    ∫d^4q/(2  π )^4    =    4  π/(2  π )^4∫^Λ_0   dq   q^3  ∫^1_-1sinϕ_1   d ( cosϕ_1 ) 
        =    1/4  π^3∫^Λ_0   dq   q^3  ∫^1_-1 √( 1 - x^2)   d x  ⋯

where Λ stands for the cutoff introduced to regulate the theory.






§ 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 <cit.>,
we provide the fits used in both works with the gluon propagator, the ghost propagator and the quark wave function curves renormalised
at μ = 4.3 GeV within the MOM scheme. On Fig. <ref> the curves referred to as JHEP are those of <cit.>,
while those designated as NEW are the curves mentioned in Secs. <ref> and <ref>.
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.









99

Marciano:1977su
  W. J. Marciano and H. Pagels,
  
  Phys. Rept.  36 (1978) 137.
  doi:10.1016/0370-1573(78)90208-9
  
Alkofer:2000wg
  R. Alkofer and L. von Smekal,
  
  Phys. Rept.  353 (2001) 281
  doi:10.1016/S0370-1573(01)00010-2
  [hep-ph/0007355].
  
Fischer:2006ub
  C. S. Fischer,
  
  J. Phys. G 32 (2006) R253
  doi:10.1088/0954-3899/32/8/R02
  [hep-ph/0605173].
  
Binosi:2009qm
  D. Binosi and J. Papavassiliou,
  
  Phys. Rept.  479 (2009) 1
  doi:10.1016/j.physrep.2009.05.001
  [arXiv:0909.2536 [hep-ph]].  
  
Bermudez:2017bpx
  R. Bermudez, L. Albino, L. X. Gutiérrez-Guerrero, M. E. Tejeda-Yeomans and A. Bashir,
  
  Phys. Rev. D 95 (2017) no.3,  034041
  doi:10.1103/PhysRevD.95.034041
  [arXiv:1702.04437 [hep-ph]].

Ball:1980ay
  J. S. Ball and T. W. Chiu,
  
  Phys. Rev. D 22 (1980) 2542.
  doi:10.1103/PhysRevD.22.2542

Curtis:1990zs
  D. C. Curtis and M. R. Pennington,
  
  Phys. Rev. D 42 (1990) 4165.
  doi:10.1103/PhysRevD.42.4165

Bashir:2011dp
  A. Bashir, R. Bermudez, L. Chang and C. D. Roberts,
  
  Phys. Rev. C 85 (2012) 045205
  doi:10.1103/PhysRevC.85.045205
  [arXiv:1112.4847 [nucl-th]].  

Qin:2013mta
  S. X. Qin, L. Chang, Y. X. Liu, C. D. Roberts and S. M. Schmidt,
  
  Phys. Lett. B 722 (2013) 384
  doi:10.1016/j.physletb.2013.04.034
  [arXiv:1302.3276 [nucl-th]].
  
Aslam:2015nia
  M. J. Aslam, A. Bashir and L. X. Gutierrez-Guerrero,
  
  Phys. Rev. D 93 (2016) no.7,  076001
  doi:10.1103/PhysRevD.93.076001
  [arXiv:1505.02645 [hep-th]].

Binosi:2016wcx
  D. Binosi, L. Chang, J. Papavassiliou, S. X. Qin and C. D. Roberts,
  
  Phys. Rev. D 95 (2017) no.3,  031501
  doi:10.1103/PhysRevD.95.031501
  [arXiv:1609.02568 [nucl-th]].
    
Alkofer:2008tt
  R. Alkofer, C. S. Fischer, F. J. Llanes-Estrada and K. Schwenzer,
  
  Annals Phys.  324 (2009) 106
  doi:10.1016/j.aop.2008.07.001
  [arXiv:0804.3042 [hep-ph]].

Aguilar:2010cn
  A. C. Aguilar and J. Papavassiliou,
  
  Phys. Rev. D 83, 014013 (2011).
  
  

Aguilar:2011xe
  A. C. Aguilar, D. Ibañez, V. Mathieu and J. Papavassiliou,
  
  Phys. Rev. D 85 (2012) 014018
  doi:10.1103/PhysRevD.85.014018
  [arXiv:1110.2633 [hep-ph]].
  
Windisch:2012de
  A. Windisch, M. Hopfer and R. Alkofer,
  
  Acta Phys. Polon. Supp.  6 (2013) 347
  doi:10.5506/APhysPolBSupp.6.347
  [arXiv:1210.8428 [hep-ph]].
  
Hopfer:2013via
  M. Hopfer and R. Alkofer,
  
  Acta Phys. Polon. Supp.  6 (2013) no.3,  929
  doi:10.5506/APhysPolBSupp.6.929
  [arXiv:1304.4360 [hep-ph]].
  
Aguilar:2014lha
  A. C. Aguilar, D. Binosi, D. Ibañez and J. Papavassiliou,
  
  Phys. Rev. D 90 (2014) no.6,  065027
  doi:10.1103/PhysRevD.90.065027
  [arXiv:1405.3506 [hep-ph]].
  
Pelaez:2015tba
  M. Peláez, M. Tissier and N. Wschebor,
  
  Phys. Rev. D 92 (2015) no.4,  045012
  doi:10.1103/PhysRevD.92.045012
  [arXiv:1504.05157 [hep-th]].    

Binosi:2014aea
  D. Binosi, L. Chang, J. Papavassiliou and C. D. Roberts,
  
  Phys. Lett. B 742 (2015) 183
  doi:10.1016/j.physletb.2015.01.031
  [arXiv:1412.4782 [nucl-th]].

Sanchis-Alepuz:2015qra
  H. Sanchis-Alepuz and R. Williams,
  
  Phys. Lett. B 749 (2015) 592
  doi:10.1016/j.physletb.2015.08.067
  [arXiv:1504.07776 [hep-ph]].    
  
Williams:2015cvx 
  R. Williams, C. S. Fischer and W. Heupel,
  
  Phys. Rev. D 93 (2016) no.3,  034026
  doi:10.1103/PhysRevD.93.034026
  [arXiv:1512.00455 [hep-ph]].

Aguilar:2016lbe
  A. C. Aguilar, J. C. Cardona, M. N. Ferreira and J. Papavassiliou,
  
  Phys. Rev. D 96 (2017) no.1,  014029
  doi:10.1103/PhysRevD.96.014029
  [arXiv:1610.06158 [hep-ph]].
  
Aguilar:2018epe
  A. C. Aguilar, J. C. Cardona, M. N. Ferreira and J. Papavassiliou,
  
  Phys. Rev. D 98 (2018) no.1,  014002
  doi:10.1103/PhysRevD.98.014002
  [arXiv:1804.04229 [hep-ph]].  

Skullerud:2001aw
  J. Skullerud, D. B. Leinweber and A. G. Williams,
  
  Phys. Rev. D 64 (2001) 074508
  doi:10.1103/PhysRevD.64.074508
  [hep-lat/0102013].  

Skullerud:2003qu
  J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber and A. G. Williams,
  
  JHEP 0304 (2003) 047
  doi:10.1088/1126-6708/2003/04/047
  [hep-ph/0303176].

Kizilersu:2006et
  A. Kizilersu, D. B. Leinweber, J. I. Skullerud and A. G. Williams,
  
  Eur. Phys. J. C 50 (2007) 871
  doi:10.1140/epjc/s10052-007-0250-6
  [hep-lat/0610078].  

Oliveira:2016muq
O. Oliveira, A. Kızılersu, P. J. Silva, J.-I. Skullerud, A. Sternbeck, A. G. Williams,
Acta Phys. Polon. Supp. 9, 363 (2016) [arXiv:1605.09632]

Rojas:2013tza
  E. Rojas, J. P. B. C. de Melo, B. El-Bennich, O. Oliveira and T. Frederico,
  
  JHEP 1310 (2013) 193
  doi:10.1007/JHEP10(2013)193
  [arXiv:1306.3022 [hep-ph]].
  
Oliveira2018a
  O. Oliveira, T. Frederico, W. de Paula and J. P. B. C. de Melo,
  
  Eur. Phys. J. C 78 (2018) no.7,  553
  doi:10.1140/epjc/s10052-018-6037-0
  [arXiv:1807.00675 [hep-ph]].

Chang:2011ei
  L. Chang and C. D. Roberts,
  
  Phys. Rev. C 85 (2012) 052201
  doi:10.1103/PhysRevC.85.052201
  [arXiv:1104.4821 [nucl-th]].

Bashir:1994az 
  A. Bashir and M. R. Pennington,
  
  Phys. Rev. D 50, 7679 (1994)
  doi:10.1103/PhysRevD.50.7679  
  
Maris99a
  P. Maris and P. C. Tandy,
  
  Phys. Rev. C 60 (1999) 055214
  doi:10.1103/PhysRevC.60.055214
  [nucl-th/9905056].


Bashir:2012fs
  A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y. X. Liu, C. D. Roberts and P. C. Tandy,
  
  Commun. Theor. Phys.  58 (2012) 79
  doi:10.1088/0253-6102/58/1/16
  [arXiv:1201.3366 [nucl-th]].

Mojica:2017tvh
  F. F. Mojica, C. E. Vera, E. Rojas and B. El-Bennich,
  
  Phys. Rev. D 96 (2017) no.1,  014012
  doi:10.1103/PhysRevD.96.014012
  [arXiv:1704.08593 [hep-ph]].


Serna:2018dwk
  F. E. Serna, C. Chen and B. El-Bennich,
  
  arXiv:1812.01096 [hep-ph].


Qin:2011dd
  S. x. Qin, L. Chang, Y. x. Liu, C. D. Roberts and D. J. Wilson,
  
  Phys. Rev. C 84 (2011) 042202
  doi:10.1103/PhysRevC.84.042202
  [arXiv:1108.0603 [nucl-th]].
  
Kizilersu:1995iz
  A. Kizilersu, M. Reenders and M. R. Pennington,
  
  Phys. Rev. D 52 (1995) 1242
  doi:10.1103/PhysRevD.52.1242
  [hep-ph/9503238].
      
Davydychev A.I.Davydychev, A. Osland and L. Saks. Phys. Rev. D 63, 014022 (2000).

  
 Slavnov:1972fg 
  A. A. Slavnov,
  
  Theor. Math. Phys.  10, 99 (1972)
  [Teor. Mat. Fiz.  10, 153 (1972)].
  
   
 Taylor:1971ff 
  J. C. Taylor,
  
  Nucl. Phys. B 33, 436 (1971).
  

Oliveira:2018lln
  O. Oliveira, P. J. Silva, J. I. Skullerud and A. Sternbeck,
  
  arXiv:1809.02541 [hep-lat].
  
Skullerud:2000un 
  J. I. Skullerud and A. G. Williams,
  
  Phys. Rev. D 63, 054508 (2001)
  doi:10.1103/PhysRevD.63.054508
  [hep-lat/0007028].

Dudal:2010tf
D. Dudal, O. Oliveira, N. Vandersickel, Phys. Rev. D 81, 074505 (2010) [arXiv:1002.2374]

Ayala:2012pb
A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, J. Rodríguez-Quintero, Phys. Rev. D 86, 074512 (2012) [arXiv:1208.0795]

Duarte:2016iko
A. G. Duarte, O. Oliveira, P. J. Silva, Phys. Rev. D 94, 014502 (2016) [arXiv:1605.00594]

Dudal:2018cli
  D. Dudal, O. Oliveira and P. J. Silva,
  
  Annals Phys.  397 (2018) 351
  doi:10.1016/j.aop.2018.08.019
  [arXiv:1803.02281 [hep-lat]].

Dudal:2012zx 
  D. Dudal, O. Oliveira and J. Rodriguez-Quintero,
  
  Phys. Rev. D 86, 105005 (2012)
  doi:10.1103/PhysRevD.86.105005, 10.1103/PhysRevD.86.109902
  [arXiv:1207.5118 [hep-ph]].
  
Aguilar:2018csq
  A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo and J. Papavassiliou,
  
  arXiv:1811.08961 [hep-ph].  
  
Tikhonov95
A.N. Tikhonov, A.V. Goncharsky, V.V. Stepanov, A.G. Yagola, Numerical Methods for the Solution of Ill Posed Problems, Kluwer Academic Publishers, 1995    

Eichmann:2016yit
  G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer and C. S. Fischer,
  
  Prog. Part. Nucl. Phys.  91 (2016) 1
  doi:10.1016/j.ppnp.2016.07.001
  [arXiv:1606.09602 [hep-ph]].
      




