Nonperturbative study of the four gluon vertex

In this paper we study the nonperturbative structure of the SU(3) four-gluon vertex in the Landau gauge, concentrating on contributions quadratic in the metric. We employ an approximation scheme where"one-loop"diagrams are computed using fully dressed gluon and ghost propagators, and tree-level vertices. When a suitable kinematical configuration depending on a single momentum scale $p$ is chosen, only two structures emerge: the tree-level four-gluon vertex, and a tensor orthogonal to it. A detailed numerical analysis reveals that the form factor associated with this latter tensor displays a change of sign (zero-crossing) in the deep infrared, and finally diverges logarithmically. The origin of this characteristic behavior is proven to be entirely due to the masslessness of the ghost propagators forming the corresponding ghost-loop diagram, in close analogy to a similar effect established for the three-gluon vertex. However, in the case at hand, and under the approximations employed, this particular divergence does not affect the form factor proportional to the tree-level tensor, which remains finite in the entire range of momenta, and deviates moderately from its naive tree-level value. It turns out that the kinematic configuration chosen is ideal for carrying out lattice simulations, because it eliminates from the connected Green's function all one-particle reducible contributions, projecting out the genuine one-particle irreducible vertex. Motivated by this possibility, we discuss in detail how a hypothetical lattice measurement of this quantity would compare to the results presented here, and the potential interference from an additional tensorial structure, allowed by Bose symmetry, but not encountered within our scheme.


I. INTRODUCTION
Of all elementary vertices that appear in the QCD Lagrangian, the four-gluon vertex is the most poorly understood. From the point of view of continuum studies, this fact may be regarded as a consequence of the enormous proliferation of allowed tensorial structures, generated by the presence of four color and four Lorentz indices. This difficulty, in turn, complicates considerably the extraction of reliable nonperturbative information from the corresponding Schwinger-Dyson equation (SDE). In addition, even gauge-technique inspired Ansätze [1][2][3][4] are extremely difficult to implement, due to the complicated structure of the Slavnov-Taylor identity that this vertex satisfies in the linear covariant (R ξ ) gauges (see, e.g. [5]). Thus, the analytic studies dedicated to this vertex are very scarce, furnishing information only at the level of one-loop perturbation theory [6,7], or involving generic constructions in the context of the pinch technique [8], or privileged quantization schemes, such as the background field method [9,10].
From the point of view of lattice simulations, the situation is simpler, in the sense that, to the best of our knowledge, no simulations of the four-gluon vertex have been performed, for any kinematic configuration. This is to be contrasted with the corresponding status of all other vertices, namely the quark-gluon, the ghost-gluon, and three-gluon vertex, which have been studied on the lattice, at least for some special choices of their momenta [11][12][13][14][15][16].
In the present work, we carry out a preliminary nonperturbative study of the one-particle irreducible (1-PI) part of the four-gluon vertex, denoted by Γ abcd µνρσ , motivated by recent developments in our understanding of the QCD nonperturbative dynamics of the two-and three-point sectors in the Landau gauge. Specifically, a precise nonperturbative connection between the masslessness of the ghost, the detailed shape of the gluon propagator in the deep infrared (IR), and the IR divergences observed in certain kinematic limits of the threegluon vertex, has been put forth in [17] (see also [18,19] for related contributions). This detailed study led to the conjecture that any purely gluonic n-point function will display the same kind of behavior, given that ghost loops 1 appear in all of them (and, hence, the associated IR logarithmic divergence in d = 4). Clearly, the confirmation of this expectation 1 We refer to ghost loops that exist already at the one-loop level. Ghost loops nested within gluon loops do not produce this particular effect, because the additional integrations over virtual momenta soften the IR divergence.
at the level of the four-gluon vertex would put our understanding of this specific IR effect on rather solid ground. In particular, it would be important to establish, even within an approximate scheme, the type of tensorial structures that will be associated with this particular divergence.
In order to simplify the calculation as much as possible without compromising its main objective, we have chosen a particularly simple configuration of the external momenta, in which a single momentum scale (p) appears, and the flow in the four legs is chosen to be (p, p, p, −3p); this has the advantage of giving rise to loop integrals that are symmetric under the crossing of external legs thus reducing the amount of diagrams one needs to evaluate. We hasten to emphasize that the aforementioned momentum configuration has been first considered in [20], in the context of the so-called "scaling" solutions [21]. Instead, our analysis will be carried out using an IR finite gluon propagator ∆ and ghost dressing function F , in conformity with the results obtained from a plethora of large-volume lattice simulations [22][23][24][25][26][27][28], as well as a variety of analytic approaches [21,. Specifically, we will consider a simplified version of the so-called "one-loop dressed" approximation, where one computes the one-loop diagrams with fully dressed gluon and ghost propagators, but with tree-level (undressed) vertices (the case with dressed ghost-vertices only is also presented).
Notice that this approach, although SDE-inspired, differs significantly from a typical SDE study, mainly because it does not involve the solution of an integral equation for the unknown form factors; instead, the form factors are simply extracted from the dressed diagrams mentioned above. In that sense, it may be thought of as a "lowest order" SDE approximation, where one simply substitutes tree-level values for all vertex form factors appearing inside diagrams. This particular method (and variations thereof) has been employed in the context of other vertices, furnishing results that compare favorably with the lattice [17,51,52]; of course, its effectiveness can only be justified a-posteriori (i.e., comparing with the lattice), given that there is no rigorous way of estimating the errors introduced by the omitted terms.
If one concentrates on the nonperturbative contributions that are quadratic in the metric, in the case of SU(3) only two independent tensorial structures emerge: the one associated with the tree-level four-gluon vertex (indicated by Γ abcd(0) µνρσ ), and a second one (denoted with G abcd µνρσ ) which is totally symmetric in both Lorentz as well color indices (and therefore orthogonal in both spaces to the tree-level term). It turns out that the aforementioned divergences are entirely proportional to this latter tensorial structure, with no contribution to the treelevel tensor Γ (0) . Therefore, one finds that within the one-loop dressed approximation we employ, G will carry all the IR divergences, whilst Γ (0) contains all the ultraviolet (UV) divergences, as required by the renormalizability of the theory. These findings clearly deviate from the patterns observed in the case of the three-gluon vertex, where the form factors proportional to the tree-level vertex, in addition to containing the UV divergences, were also affected by this particular IR divergence (displaying the associated "zero crossing").
In addition, the deviation of the form factor associated to the tree-level tensor Γ (0) from 1, namely its tree-level value, is relatively modest. In particular, when the ingredients used in its calculation are renormalized at µ = 4.3 GeV, its highest value, located at about 500 MeV, is 1.5.
The results obtained are further discussed in the specialized context of a possible future lattice simulation of the connected part of this vertex, to be denoted by C abcd µνρσ . It turns out that the momentum configuration (p, p, p, −3p) eliminates all contributions to C from oneparticle reducible (1-PR) graphs, thus isolating only Γ abcd µνρσ , without any "contamination" from lower-order Green's functions. In addition, an analysis based on Bose symmetry arguments alone, reveals that a third tensor structure, denoted by X ′abcd µνρσ , is in principle allowed; evidently, the form factor associated with this tensor vanishes within the one-loop dressed approximation that we employ. It is likely, however, that this particular property will not persist in a complete nonperturbative computation, as the one provided by lattice simulations. Therefore, under the assumption that such a structure might eventually emerge, we describe how to express the complete set of form factors characterizing Γ abcd µνρσ in terms of the standard lattice ratios R, used in studies of the three-gluon vertex [15,16].
The article is organized as follows. In Sect. II we introduce our notation, review the relevant tensor decomposition, and recall some identities particular to the SU(3) gauge group. Next, in Sect. III we carry out the calculation of the one-loop dressed diagrams in the simplified setting where all the external momenta are set to zero. This will prove to be a very useful exercise, as it will allow to determine the tensorial structures that appear, and in particular establish that the divergent part coming from ghost loops is entirely proportional to the G abcd µνρσ tensor alone. Then, in Sect. IV we carry out the calculation in the (p, p, p, −3p) momentum configuration. After manipulating all diagrams analytically (Sect. IV A), we evaluate numerically all the contributions obtained, using (quenched) lattice results as input for the gluon and ghost two-point sectors (Sect. IV B). Finally, in Sect. IV C we show how our results can be related to quantities customarily studied on the lattice. Specifically, we prove that the special momentum configuration chosen for our study has the property of isolating the 1-PI contribution to the connected four-gluon Green's function. Then, assuming the most general tensor decomposition of this vertex in terms of tensors allowed by Bose symmetry, we show what would be the best choice of the ratios R. The paper ends with Sect. V, where we draw our conclusions, and two Appendices. In the first, we carry out a general analysis of the tensor structures (quadratic in the metric) that are allowed by Bose symmetry, paying particular attention to the case (p, p, p, −3p). Finally, Appendix B collects some lengthy expressions appearing in our analytical calculations.

II. GENERALITIES ON THE FOUR-GLUON VERTEX
The 1-PI four-gluon vertex will be denoted by the expression (all momenta entering) At tree-level one has Γ abcd(0) µνρσ = f adr f cbr (g µρ g νσ − g µν g ρσ ) + f abr f rdc (g µσ g νρ − g µρ g νσ ) + f acr f dbr (g µσ g νρ − g µν g ρσ ), (2.2) where f abc are the real and totally antisymmetric SU(N) structure constants, satisfying the normalization condition so that the generators of the adjoint representation are given by In Fig. 1 we show the conventions of momenta and Lorentz/color indices used throughout this paper.
It is elementary to verify the validity of this symmetry for the tree-level vertex Γ abcd(0) µνρσ . It is clear that the fully dressed Γ abcd µνρσ is characterized, in general, by a vast proliferation of the tensorial structures (138 for general kinematics [7]); of course, as we will see, Bose symmetry imposes restrictions on the structure of the possible form factors composing Γ abcd µνρσ (p 1 , p 2 , p 3 , p 4 ). At the level of the rank-4 Minkowski tensors, the structures allowed are terms quadratic in the metric, linear in the metric and quadratic in the momenta, and quartic in momenta; schematically one has then the structures gg; gpq; pqrs. (2.5) At the level of the rank-4 color tensors the situation is considerably more complex, since, in addition to terms quadratic in f or δ, the real and totally symmetric tensors d abc will also emerge. Thus, in principle one has 15 allowed structures of the schematic type However, these tensors are related by a set of 6 identities [6], namely f abr d cdr + f acr d dbr + f adr d bcr = 0, (2.8) and two independent permutation for each, a fact that reduces the number of required tensors down to 9.
Of course, due to practical limitations, one must restrict the present study to a considerably more reduced (but physically relevant) subset of the full Lorentz and color tensorial basis mentioned above. Specifically, as was done in [6], we only consider terms quadratic in the metric tensor g µν , namely terms proportional to g µν g ρσ , g µρ g νσ , and g µσ g νρ , neglecting terms quadratic and quartic in the momenta. Thus, a priori, for a general SU(N) gauge group, one has 9 × 3 = 27 possible combinations. Furthermore, we will directly specialize our analysis to the case N = 3, where the additional identity can be used, thus reducing the number of tensorial combinations down to 24.
However, it turns out that, within the one-loop dressed approximation and the kinematical configuration that we will employ (see Fig. 2 for the 18 diagrams appearing in this case), the color tensors reduce finally to the two structures appearing in the conventional one-loop calculation of this vertex (for N = 3), namely the tree-level tensor Γ (0) defined in Eq. (2.2), and the totally symmetric (both in Minkowski and color space) tensor G abcd µνρσ = (δ ab δ cd + δ ac δ bd + δ ad δ bc ) (g µν g ρσ + g µρ g νσ + g µσ g νρ )

Rµνρσ
. (2.10) In particular, notice that since the two tensors are orthogonal in both spaces the prefactors multiplying them can be unambiguously identified 2 .
Let us finally point out that, in SU(3), one has the additional useful formula Γ abcd(0) µνρσ + G abcd µνρσ = 2X abcd µνρσ , (2.12) where we have defined the combination X abcd µνρσ = δ ab δ cd + 3 2 d abr d cdr g µν g ρσ + δ ac δ bd + 3 2 d acr d bdr g µρ g νσ Our analysis of the four-gluon vertex will be carried out in the Landau gauge, where the study of the lower Green's functions (such as gluon and ghost propagator, ghost-gluon vertex and three-gluon vertex) has been traditionally carried out, both in the continuum as well as on the lattice. In this particular gauge the full gluon propagator takes the form i∆ µν (q) = −iP µν (q)∆(q 2 ); P µν (q) = g µν − q µ q ν /q 2 , (2.14) while the ghost propagator, D(q 2 ), and its dressing function, F (q 2 ), are related by Evidently, both ∆(q 2 ) and D(q 2 ) constitute crucial ingredients for the calculations of the four-gluon vertex that follows. It is therefore useful to briefly review some of their IR features that are most relevant to the present work. Specifically, both large-volume lattice simulations and a plethora of continuous nonperturbative studies, carried out both in SU(2) and in SU(3), converge to the conclusion that the function ∆(q 2 ) reaches a finite (nonvanishing) value in the IR. Moreover, the nonperturbative ghost propagator remains "massless", and displays no IR enhancement, since its dressing function F (q 2 ) saturates in the deep IR to a finite value. As we will see in what follows, the aforementioned features have far reaching consequences for the IR behavior of the four-gluon vertex. Specifically, as happens with the tree-gluon vertex, the masslessness of the ghost-loops contributing to Γ abcd µνρσ produces a logarithmic IR divergence. What is, however, qualitatively distinct compared to the three-gluon case, is that, at least within the approximation scheme that we employ, this particular divergence does not manifest itself in the part proportional to Γ (0) , but rather in the orthogonal combination G.

III. VANISHING EXTERNAL MOMENTA
In this section we consider the simplest possible kinematic case, where all the momenta of the external gluons are set to zero (p 1 = p 2 = p 3 = p 4 = 0).

A. The calculation
Since we do no consider the contribution of quark-loops (pure Yang-Mills theory), the only representation that appears in our problem is the adjoint, whose explicit realization is given in Eq. (2.4).
For the various integrals appearing in this calculation we will employ the standard text- book results where R µνρσ has been defined in Eq. (2.10), and the integral measure is k = µ ǫ d d k/(2π) d , with d = 4 − ǫ the space-time dimension 3 and µ the 't Hooft mass.
There are two particular tensorial structures that appear in a natural way in the calculations of the graphs shown in Fig. 2, namely Then, using the relation 4

3)
3 Notice that we set d = 4 − ǫ instead of d = 4 + 2ǫ used in [6]. 4 Note also the particular property Tr adj ( and G, Turning to the explicit calculation of the one-loop dressed diagrams of Fig. 2, the (six) ghost boxes give the result which, with the aid of the formulas (3.4) introduced above, may be written in the simple Since the ghost dressing function F is known to saturate in the IR, the integral above diverges logarithmically in the IR; however Eq. (3.6) shows that this divergence does not contribute to the structures proportional to the tree-level tensor Γ (0) . Even though this result has been derived in a simplified setting, it will persist within the one-loop dressed approximation employed here. Therefore, we arrive at the important conclusion that the IR divergent terms originating from the ghost loops would be completely missed, if one were to consider only the form factor proportional to the tree-level tensor Γ (0) .
We next consider the (three) gluon boxes; as the adjoint traces will be the same as those appearing in the ghost case above, we obtain that also the one-loop dressed gluon boxes do not contribute to the tree-level tensor structure. In particular, we get Notice that, unlike the case of the ghost boxes treated above, the integral appearing in Eq. (3.7) is convergent in the IR, because the gluon propagator reaches a finite value in that limit.
We now turn to the (six) triangle diagrams. After some straightforward algebraic manip- which, after using the identities (3.4), can be cast in the form Finally, we are left with the (three) fishnet diagrams. One finds, similarly to what happens with the triangle diagrams, The identities (3.4) allow us to express the result in its final form, namely The results obtained are conveniently summarized in Table I.

B. Perturbative analysis
At this point one may explore the qualitative behavior of the two contributions obtained above within a setting inspired by one-loop perturbation theory, but supplemented by a set of mass scales, which prevent the resulting expressions from diverging in the IR. Specifically, if one were to simply set F (k 2 ) and ∆(k 2 ) to their strict perturbative values (1 and 1/k 2 , respectively) the four integrals appearing in the second column of Table I reduce to a single integral, namely k 1 k 4 . At this point, it is easy to verify that, when d = 4, the total contribution proportional to G abcd µνρσ vanishes, given that the sum of the coefficients appearing on the fourth column adds up to zero.
However, given that the integral k 1 k 4 is both IR and UV divergent, it is preferable to introduce a distinction between the two type of divergences. To accomplish this, we proceed as follows. Given that the (Euclidean) gluon propagator (in the Landau gauge) is known to be finite in the IR (a feature that can be self-consistently explained through the dynamical generation of an effective gluon mass), for the purposes of this simple calculation one may approximate ∆(k 2 ) by 1/(k 2 + m 2 ). This replacement makes the integrals k k 4 ∆ 4 (k 2 ), Table I IR finite; of course, they still diverge logarithmically in the UV. Regarding the integral k F 4 (k 2 ) k 4 , it is known that the ghost remains nonperturbatively massless, a fact that leads to a genuine IR divergence; in order to control it, we will introduce an artificial mass scale, denoted by λ 2 . Thus, the integral corresponding to k Let us emphasize at this point that even though at the formal level both m 2 and λ 2 serve as IR regulators, there is a profound physical difference between the two: m 2 constitutes a simplified realization of a true physical phenomenon, namely the IR saturation of the gluon propagator, while λ 2 is an artificial scale, introduced as a regulator of a quantity (the ghost propagator) that is genuinely massless. Consequently, in order to recover the physically relevant (albeit simplified) limits, m 2 will be kept at some fixed nonvanishing value, while λ 2 will be sent to zero.
The above considerations motivate the introduction of a particular integral, namely where µ is the 't Hooft mass, and γ the Euler-Mascheroni constant. Evidently, depending on the case that one considers, M 2 = m 2 or M 2 = λ 2 .
In particular, after the replacements mentioned above, the integrals in Table I can be expressed in terms of I(M 2 ) as follows  Table I and obtain, within this perturbative scheme, the coefficients multiplying Γ (0) and G, to be denoted by V which, after the inclusion of the tree-level term, and use of Eq. (3.14), becomes and V (1) Evidently, all dependence on 1/ǫ is contained in the coefficient multiplying Γ (0) , while the coefficient of G is completely free of such terms, exactly as one would expect from the renormalizability of the theory. Indeed, given that the term G does not appear in the original Lagrangian, a divergence of this type could not be renormalized away. Instead, the divergence proportional to Γ (0) will be reabsorbed in the standard way, namely through the introduction of the appropriate vertex renormalization constant, to be denoted by Z 4 .
Specifically, one obtains the renormalized vertex Γ R from its unrenormalized counterpart Γ 0 through the condition (suppressing all indices) Of course, the exact form of the Z 4 and the resulting Γ R depend on the renormalization scheme chosen. In particular, in the minimal subtraction (MS) scheme one would simply which, upon multiplication with the V Γ (0) (0) of Eq. (3.17) will give (keeping up to terms of coincides with the part proportional to 1/ǫ of the corresponding expression given in (3.10) of [6] (in the Landau gauge, and for N = 3).
If one were instead to renormalize in the minimal subtraction (MOM) scheme, as is customary in lattice simulations and SDE studies, one would need to introduce a renormalization point, µ R , and demand that at that point the value of the renormalized vertex reduces to its tree-level value. For instance, as in [6], the completely symmetric choice such that (schematically) Of course, for the case at hand, since the vertex has been computed only for vanishing momenta, one cannot implement a MOM-type scheme. However, in order to get a sense of the general trend that one might expect from a general calculation, we may assume that the subtraction point lies sufficiently far in the UV. Then, for a representative large Euclidean momentum P , the qualitative behaviour of the form factor may be approximated by and therefore, the value of V Γ (0) (0) gets renormalized to As happens typically, in the finite result the 't Hooft scale µ has been replaced by the renormalization scale µ R .
It is obvious at this point, that the above approximations require that µ 2 R > m 2 , and, consequently, since the logarithm becomes negative, V To obtain a quantitative notion of the effect, we will use lattice-inspired values for m 2 and µ 2 R ; specifically, if we identify the saturation point of the gluon propagator on the lattice with 1/m 2 , we know that, for µ R = 4.3 GeV we have that m = 375 MeV. Then, using that, for this particular Turning to the V G (0) in Eq. (3.18), we notice that, when the artificial IR cutoff λ is taken to zero, while the physical gluon mass is kept at a nonvanishing value, the logarithm diverges to +∞. Again, this coincides with the behavior found in the more complete calculation of the next section. Of course, the slope of the logarithm found in Eq. (3.18) is numerically rather suppressed when compared to the result found in the next section; however, this is to be expected, given that the function F (k 2 ), which in Eq. (3.6) is raised to the fourth power, is considerably different from 1 in the IR and intermediate momenta. Even within the one-loop dressed approximation we are employing, the calculation of the four-gluon vertex for a generic external momenta configuration (such as the one depicted in Fig. 1) is still a complex task. In addition, it is not the most expeditious way to obtain information about the IR dynamics of this vertex that could be easily contrasted with lattice simulations.
Thus, we will study a relatively simple kinematic configuration, which is obtained choosing a single momentum scale p and identifying the momentum flow (see Fig. 1) with p 1 = p 2 = p 3 = p (and hence p 4 = −3p). This kinematic configuration gives rise to loop integrals that are fully symmetric under the crossing of external legs; therefore, the crossed diagrams may be obtained from the original ones through simple permutations of the color and Lorentz indices.
As before, we will only consider terms that are quadratic in the metric gg. This choice, in addition to simplifying the algebraic structures considerably, corresponds precisely to the contributions that would survive on the lattice, if one were to consider the standard quantities employed in the simulations of vertices [15,16] (we will return to this point in Sect. IV C).

A. Analytical results
Consider the contribution of the ghost boxes. The aforementioned crossing property implies that the six different diagrams are proportional to the same integral. As a result, one obtains, similarly to what happens in the p = 0 case, where now It can be easily checked that as p → 0, A(p 2 ) above reduces to the A(0) of Eq. (3.6); therefore we expect that the G abcd µνρσ form factor will develop a (logarithmic) divergence in the deep IR. Next, we consider the gluon boxes. The uncrossed diagram shown in Fig. 2, yields the general expression where the integrals I i (p 2 ) are not needed for the moment. Crossed diagrams are then obtained from the above expression through the replacements (µνρσ) → (µρνσ), (abcd) → (acbd) and (µνρσ) → (νµρσ), (abcd) → (bacd). In addition, it turns out that the integrals I 1 and I 3 are equal 6 upon the momentum shifting k + 3p → −k, so that I µνρσ can be cast 6 Notice that without this equality the gluon box contributions would lie outside the subset of all possible color and Lorentz tensor structures spanned by Γ abcd(0) µνρσ and G abcd µνρσ . Moreover, observe that the realization of this equality requires shifts of the integration variable of the type k → k + p; of course, since only terms quadratic in the metric are kept, one consistently drops in the numerators terms produced by these shifts that are proportional to p and carry Lorentz indices of the external legs. in the form Thus, adding the three diagrams, one obtains 6) and the functions f i (k, p) are reported in Eq. (B1).
We next consider the triangle diagrams. In this case the six graphs can be divided in two separate classes (see Fig. 3), proportional to two independent momentum integrals, namely where J µνρσ (p 2 ) = J 1 (p 2 )(g µν g ρσ − g µρ g νσ ), K µνρσ (p 2 ) = K 1 (p 2 )g µσ g νρ + K 2 (p 2 )g µρ g νσ + K 3 (p 2 )g µν g ρσ + K 4 (p 2 )R µνρσ , where again J i (p 2 ) and K i (p 2 ) are integrals whose explicit expression is not needed at this point.

FIG. 4: (color online). The SU(3) gluon propagator (left) and ghost dressing function (right)
evaluated on the lattice [25] and the corresponding physically motivated fits we use [53]. In the case of the gluon propagator the dashed curve shows a fit featuring an inflection point the origin of which is linked to the presence of ghost loops [17]. All functions are renormalized at µ = 4.3 GeV.
where now where the h i (k, p) functions are given in Eq. (B4).
At this point, using the above results, and taking into account the definition (2.1), one has that the four-gluon vertex can be cast in the form where the "1" in V Γ (0) (p 2 ) represents the tree-level contribution.

B. Numerical results
In order to study numerically the various one-loop dressed contributions to the four-gluon vertex, let us first pass to Euclidean space by defining k 0 → ik E 4 and k j → −k E j , from which the replacement rules d 4 k → id 4 k E , k ·q → −k E ·q E and k 2 → −k 2 E follow. Next, we introduce spherical coordinates, setting x = p 2 ; y = k 2 ; z n = (k + np) 2 = n 2 x + y + 2n √ xy cos θ; At this point, all the integrals derived in our analytical calculation may be evaluated by standard integration techniques, provided that we supply as input the gluon propagator ∆ and the ghost dressing function F .
To this end, we use physically motivated fits to the lattice data of [54], whose explicit functional form can be found in [53]. The agreement of these fits with the corresponding lattice data at the renormalization scale µ = 4.3 GeV is shown in Fig. 4. For the case of the gluon propagator we also show a fit displaying the inflection point that must appear due to the presence of divergent ghost loops [17]; the results obtained are practically independent from the implementation of this feature in the gluon propagator.
In Fig. 5 and Fig. 6 we plot, respectively, the contributions of the various diagrams to V Γ (0) and V G , together with their total sum (in the former case, all terms have been subtractively renormalized within the MOM scheme, at µ = 4.3 GeV, in accordance with Eq. (3.23)). As already mentioned, ghost boxes will not contribute to V Γ (0) , which is entirely made up of gluonic contributions, all of them saturating in the IR (see Fig. 5, again). The contribution of the gluon boxes is negligible; indeed, as p → 0 it vanishes, as we know it should from the zero external momenta case (see Table I). The triangle terms feature a bump of opposite sign, while the fishnet is negative. Adding everything up, one obtains the shape shown by the black line of Fig. 5. Notice that at zero momentum we obtain the value V In the case of V G the situation is completely different (Fig. 6). Gluon contributions are again saturating in the IR; however, in this case, the ghost boxes take over below few hundreds MeV 2 , driving V G to an IR logarithmic divergence. In fact, the IR behavior is perfectly described by the function a log x + b with a = −0.187 and b = −1.989. As far as the remaining diagrams are concerned, gluon boxes are negative in this case; in addition, they are almost perfectly cancelled by the two triangle contributions, which (contrary to the previous case) have now the same sign. When the negative contribution from the fishnet diagrams is finally added, one obtains the shape shown by the black line of Fig. 6.
It is important to notice that V G displays a zero crossing, a feature that is also present in the R ratio defined in the case of the three-gluon vertex [17][18][19]. The location of the  the most general tensorial structure decomposition of Γ α is given by Evidently, at tree-level, one has A (0) = 1 and B (0) = 0.
To be sure, the form factors A, B are not known for arbitrary momenta. However, in the soft gluon limit (p → 0), which is the most relevant to our purposes, Eq. (4.21) reduces to Γ α (k, 0, −k) = A(k 2 )k α ; A(k 2 ) = A(k, 0, −k), (4.22) and the form factor A(k 2 ) has been studied both in the continuum [51,55,56] and on the lattice, both for SU(2) [13] and SU(3) [14]; the results obtained in [51] and [14] are shown in the left panel of Fig. 7. As can be seen, in this limit A develops a sizeable peak around 800 MeV, approaching its tree-level value for both IR as well as UV momenta.
The idea is then to replace all ghost-gluon vertices appearing in a generic ghost box diagram by Eq.  (4.23) and we see that 1-PR diagrams constructed from lower order functions will spoil in general the possibility of isolating the genuine 1-PI contribution to C. In order to address the second problem, observe that, in the Landau gauge, the only rank-2 Minkowski tensor allowed is the transverse projector, while, in general, the allowed rank-3 tensors for the three-gluon vertex are either linear in both the metric and momenta (gp), or cubic in momenta (pqr). This means, in turn, that in the momentum configuration (p 1 , p 2 , p 3 , p 4 ) = (p, p, p, −3p), each propagator appearing in the decomposition Eq. (4.23) will have an accompanying projector P (p); as a result, all 1-PR contributions will vanish, and one is left with At this point, the scalar factors ∆ can be factored out by defining the lattice R ratio in the standard way [15,16], namely projecting C on a suitable tensor T abcd µνρσ , and normalizing the resulting expression. Specifically, one writes T abcd µνρσ P µα P νβ P ργ P σδ T abcd αβγδ = T abcd µνρσ P µα P νβ P ργ P σδ Γ abcd αβγδ gg T abcd µνρσ P µα P νβ P ργ P σδ T abcd αβγδ . (4.25) One then usually chooses T to coincide with the tree-level vertex Γ (0) , so that any deviation of R from 1 signals the onset of quantum (nonperturbative) effects.
For the case of the four gluon vertex, however, additional care is needed, depending on the particular property that one attempts to expose 7 . Indeed, the general analysis based on Bose symmetry of Appendix A, reveals that in the momentum configuration under scrutiny there are at most three possible tensor structures (proportional to gg) contributing to the full four-gluon vertex, i.e., one has 8 Thus, the complete structure of the four gluon vertex can be obtained by defining the three different ratios corresponding to setting T = Γ (0) , G and X ′ in the definition (4.25).
In addition, as explained in Appendix A, Bose symmetry alone does not unambiguously fix the tensor X ′ , whose exact form depends on the choice of the "basis" that spans this Ideally one would like to fix s in a way such that the resulting X ′ be orthogonal to both Γ (0) and G, that is by requiring This is however not possible, as X ′ can be rendered orthogonal to either Γ (0) (for s = 0) or G (for s = 1/3), but not both. In these cases one has There are at least two reasons to prefer the second choice over the first. To begin with, recall that, according to our general analysis, the origin of the divergence in the vertex form factors is clearly associated with the masslessness of the full ghost propagator; consequently, V G will continue to be divergent even in the full nonperturbative setting provided by a lattice calculation. At the same time, we expect V Γ (0) to be finite, as the diagrams contributing to it are "protected" by the effective gluon mass; this means, in turn, that in the s = 0 basis both R G and R X ′ will diverge, even in the case of a finite V X ′ . In the s = 1/3 basis, a lattice calculation would instead find that the only IR divergent ratio would be R G , immediately signalling a finite V Γ (0) and V X ′ form factors. In addition, observe that Eq. (4.28) implies 8 Obviously, within the one-loop approximation that we have employed in our calculations, one has V X ′ = 0, that the vectors Γ (0) and X ′ are very close to be orthogonal for the s = 1/3 case; for example, when projecting the full vertex along Γ (0) , the X ′ component is two orders of magnitude smaller than the Γ (0) one, and vice-versa.
A lattice measurement of the ratios R T in the s = 1/3 basis will finally yield the complete vertex form factors V T , through the formulas It turns out that ghost boxes can contribute only to the latter structure, while all the remaining diagrams (gluon boxes, triangle, and fishnet diagrams) contributes to both. Then, as massless ghost loops invariably lead to the presence of an IR divergence in the corresponding diagram [17], one expects the form factor V G (respectively, V Γ (0) ) to be IR divergent (respectively, finite). A numerical study performed, using as input the available lattice data for the gluon propagator and ghost dressing function, confirms these expectations. In addition, one finds that V G shows a zero crossing before the form factor diverges logarithmically to +∞.
It would certainly be interesting to scrutinize this issue further, and reach a definite conclusion on the way that this particular IR divergence manifests itself at the level of the four-gluon vertex. One possible direction has already been pursued here to some limited extent, namely the dressing of the ghost vertices appearing in the one-loop dressed diagrams.
One may attempt to complete this task, by also dressing the three-gluon vertices; to be sure, the tensorial structure of the three-gluon vertices is bound to lead to a proliferation of terms, which, however, may become manageable in the limit of interest, namely as all external momenta tend to zero. Even if this approach would not exhaust all the possible vertex dressing (the real problem in this context being the dressing of the four-gluon vertices appearing in the triangle diagrams), the importance of accomplishing this step would be twofold: on the one hand, one would see if the zero crossing gets pushed towards the more "favorable" momentum region p ∼ 0.1 ÷ 1 GeV, as the ghost vertex corrections seem to suggest; on the other hand, it might be possible to detect the appearance of the form factor associated to the X ′ tensor allowed by Bose symmetry, and address its IR properties (or confirm its vanishing).
Of course, the lattice could be instrumental in addressing all the aforementioned issues.
Indeed, our analysis reveals that the (p, p, p, −3p) configuration would permit the study of the 1-PI part of the connected four gluon function alone, and that the full structure of the vertex can be then reconstructed from the measurements of the standard ratios R Γ (0) , R G and R X ′ . In that sense, the main remaining difficulty to overcome is to average over a large sample of gauge configurations, in order to tame the statistical fluctuations.

Acknowledgments
The research of J. P. is supported by the Spanish MEYC under grant FPA2011-23596. for this particular separation is that the elements of different subsets do not mix when one applies the permutations dictated by the Bose symmetry of the four-gluon vertex, as we will do below.
In particular, the three basic structures that emerge when one considers only terms quadratic in the metric assume the general form and the vertex may be written as The Let us first consider the contributions to V i coming from the first subset, to be denoted by V A i . In particular, we have (V A 1 ) abcd µνρσ = g µν g ρσ [a 11 δ ab δ cd + a 12 δ ac δ bd + a 13 δ ad δ bc ], (V A 2 ) abcd µνρσ = g µρ g νσ [a 21 δ ab δ cd + a 22 δ ac δ bd + a 23 δ ad δ bc ], (V A 3 ) abcd µνρσ = g µσ g νρ [a 31 δ ab δ cd + a 32 δ ac δ bd + a 33 δ ad δ bc ].
Note that if one carries out the second obvious set of permutations, namely (c, ρ ↔ d, σ), , respectively, are automatically symmetric.
Of course, when the permutation is such that one particular V A i must transform into itself, the other two must transform one into the other. For example, when (a, µ ↔ b, ν), we Then, Bose symmetry requires that the transformed V A 2 must coincide with the original V A 3 , and vice-versa, and therefore we must have that a 22 = a 33 , and a 23 = a 31 . The repetition of this arguments leads to the conclusion that a 11 = a 22 = a 33 ≡ a, and a 12 = a 23 = a 31 ≡ a; thus, finally, after setting a ≡ a − a, we have that (V A 1 ) abcd µνρσ = g µν g ρσ [aδ ab δ cd + a(δ ab δ cd + δ ac δ bd + δ ad δ bc )], (V A 2 ) abcd µνρσ = g µρ g νσ [aδ ac δ bd + a(δ ab δ cd + δ ac δ bd + δ ad δ bc )], (V A 3 ) abcd µνρσ = g µσ g νρ [aδ ad δ bc + a(δ ab δ cd + δ ac δ bd + δ ad δ bc )].
Note that past this point, use of the remaining possible permutations imposes no further restrictions on the coefficients a and a.
A completely similar procedure may be applied to the parts of the (V i ) related to the subset {B}. In particular, one reaches the conclusion that and since the transformed (V C 1 ) must coincide with the original one, we have that c 11 = −c 11 and c 12 = −c 13 , and so where we have used the second identity of Eq. (2.8). But this last expression must remain invariant under the additional permutation (a, µ ↔ b, ν), which implies that c 13 = 0.
A this point one may specialize to the case N = 3, and use Eq. (2.9) into Eq. (A8), to write the vertex in the form with L abcd µνρσ = aE abcd µνρσ + bE ′ abcd µνρσ , E abcd µνρσ = g µν g ρσ δ ab δ cd + g µρ g νσ δ ac δ bd + g µσ g νρ δ ad δ bc , E ′ abcd µνρσ = g µν g ρσ d abr d cdr + g µρ g νσ d acr d bdr + g µσ g νρ d adr d bcr . (A13) The term L may be further manipulated, by noticing that if the condition a = 2b 3 were satisfied, then we would have that L abcd µνρσ = 2b µνρσ ] [see Eq. (2.12)]. Therefore, the most general way to rearrange this term is where we have set a = 2b 3 + c, and s represents a freely adjustable parameter that controls the weights with which the tensors E and E ′ enters the definition of the vector X ′ . Then, Γ abcd µνρσ (0, 0, 0, 0) may be decomposed in terms of the color and Lorentz vectors Γ Evidently, within the one-loop dressed approximation one has c = 0. We next turn to the case (p 1 , p 2 , p 3 , p 4 ) = (p, p, p, −3p). In this case, the general form of the V i given Eq. (A2) remains the same, but now the form factors are functions of the only available momentum scale, namely p 2 , so that a ij → a ij (p 2 ), b ij → b ij (p 2 ) and c ij → c ij (p 2 ).
In general, the presence of momenta makes the implementation of Bose symmetry more complicated, because it involves additional permutations of (p 1 , p 2 , p 3 , p 4 ). However, for the particular case at hand, the fact that the form factors can only depend on p 2 , makes these momentum permutations "inert". As a result, one arrives at exactly the same form for Γ abcd µνρσ (p, p, p, −3p) as the one given in Eq. (A12), with all coefficients converted into functions of p 2 .
Thus, after the definition a(p 2 ) ≡ a(p 2 ) − a(p 2 ), one arrives at exactly the same expression as in Eq. (A7), with the only difference that the coefficients are now functions of p 2 . The same conclusions are reached for the V B i , where b → b(p 2 ) and b → b(p 2 ), while V C i vanishes as before.

(B4)
It is straightforward but tedious to verify that, in the limit p → 0, the above expressions reduce to the corresponding results found in Sect. III.