Four Gluon Vertex from Lattice QCD

A lattice QCD calculation for the four gluon one-particle irreducible Green function in the Landau gauge is discussed. Results for some of the associated form factors are reported for kinematical configurations with a single momentum scale. Our results show that the computation of this Green function requires large statistical ensembles with 10K or larger number of gauge configurations. The simulations considered herein have a clear Monte Carlo signal for momenta up to $\sim 1$ GeV. The form factors show an hierarchy, with the form factor associated with the tree level Feynman rule being dominant and essentially constant for the range of momenta accessed. The remaining form factors seem to increase as the momentum decreases, suggesting that a possible $\log$ divergence may occur. The computed form factors are, at least, in qualitative agreement with the results obtained with continuum approaches to this vertex, when available.


I. Introduction 1
II.The Four Gluon Vertex 2 A. An orthogonal set of operators 6 B. Full Green function and lattice form factors 6 C. Lattice form factors F (i) and 1-PI form factors 7 III. Lattice setup and choice of momenta 7 IV. Lattice form factors 8 V. Summary and Conclusions

I. INTRODUCTION
The dynamical content of the QCD Green functions is of upmost importance to all hadronic physics, specially its nonperturbative contents.Indeed, among other phenomena, the interpretation of the hadron spectra, the comprehension of dynamical chiral symmetry breaking and of the quark and gluon confinement mechanisms are outside the scope of the perturbative solution of the theory.
The analysis of the QCD Green functions started looking at those functions with smaller number of external legs, namely the propagator or two point correlation functions.Their tensorial structure is simpler, as it requires smaller tensor basis, and can be guided with the help of the Slavnov-Taylor identities.All the two point QCD correlation functions have been studied in detail with non-perturbative methods and a fairly good understanding of these functions has been achieved [1][2][3][4][5][6][7].A good and coherent picture of the two point QCD correlation functions, obtained with different non-perturbative methods, has been achieved but there are issues that are still under debate.An example is their analytical structure, required in the computation of meson and baryon properties with continuum methods, that remains to be understood [8][9][10][11][12].
The next level of QCD one-particle irreducible Green functions to be studied with non-perturbative methods are those with four external legs.These include the four quark and the four gluon correlation functions.The four quark Green function is required to understand, for example, the meson spectra.The four gluon correlation function has a similar role as the four quark function but for the glueball spectra.Moreover, it also allows for the definition of a renormalization group invariant QCD charge, it contributes to the Dyson-Schwinger equation for the gluon propagator and impacts on the gluon propagator itself [59].Herein, our focus of attention is the four gluon one-particle irreducible Green function.
The non-perturbative content of the four gluon 1-PI Green function has been studied with the continuum formulation of QCD using the Dyson-Schwinger equations [60][61][62].A full description of this four leg function requires a tensor basis with a large number of elements [61,63,64], that makes the analysis of this four point correlation function difficult.Continuum studies of this function have to make simplifications to arrive at a manageable computational scenario and the calculations performed have to consider a basis with a reduced number of elements.The continuum based Dyson-Schwinger calculations identify a number of form factors that display very mild dependence with the momentum, as seems to be the case of the form factor associated with the tree level tensor structure that, nevertheless, seems to be suppressed at low momenta, and others that show log divergences in the IR region.These divergences are associated with the unprotected log's coming from the massless ghost loops in the Dyson-Schwinger equations [61].The continuum QCD methods based on the Dyson-Schwinger equation for this 1-PI Green function seem to favour such type of scenario.
The picture that emerges from the various continuum approaches to the four gluon 1-PI Green function, although being consistent at the qualitative level, also show differences that probably come from the different truncations and/or the required modelling to solve the underlying equations.The differences between the various studies provide an extra motivation to a lattice calculation of the four gluon vertex.Indeed, the interplay of the various nonperturbative methods is important to achieve a proper and sufficiently accurate description of the QCD one-particle irreducible Green function.
In this work we address the computation of the four gluon one-particle irreducible Green function within pure Yang-Mills lattice QCD simulations in the Landau gauge.Preliminary results can be found in [65].The manuscript is organized as follows.In Sec.II it is given a description of the four gluon full Green function, accessed in the lattice QCD simulations, and how it relates to the corresponding 1-PI function.Furthermore, details on the extraction of the 1-PI contribution in a simulation in the Landau gauge are also discussed there, together with the tensor basis that will be used for the calculation.The definitions of the form factors to be measured can also be found in this section.The lattice setup, the choice of the kinematic configurations to be used in the measurements and some details of the lattice calculation are given in Sec.III.The lattice form factors, including those associated with the amputated Green functions, are reported and discussed in Sec.IV.Finally in Sec.V we summarize and conclude.Some of the details of the calculation, namely a proof of the decoupling of the three gluon contribution to the four point function for a certain class of kinematical configuration and the lattice version of the continuum operators, are given in the appendices.

II. THE FOUR GLUON VERTEX
Let us start by discussing the computation of the four gluon 1-PI Green function using the continuum formulation of QCD in Minkowski spacetime.By doing so we aim to arrive at a better understanding of this complex Green function before plugging into the measurement of the Green function with lattice QCD simulations.
In the numerical Monte Carlo simulations of QCD on a lattice only the full Green functions are accessed.In real space, these Green functions are defined as the vacuum expectation values of the time order product of the fundamental fields.In QCD, the gluon N-point full Green function reads where ⟨• • • ⟩ represent the vaccuum expectation value, realised as an ensemble average over gauge configurations sampled with an appropriate action.Any other quantity, like form factors, have to be computed a posteriori by handling a set of full Green functions G       and the first step to perform is to decompose the 4-point full Green function in terms of QCD 1-PI Green functions.
For pure Yang-Mills theory the functional generator of the full Green functions Z[J] is where S[A] is the effective action associated with the Yang-Mills theory, that has contributions coming from the gauge and ghost fields, W [J] is the generating functional of the connected Green functions and J are classical external sources.It follows from the definition of the full Green function that the four gluon full Green function is given by The diagrammatic representation of ) in terms of connected Green functions is given in Fig. 1.The four gluon full Green function is a sum of a four-gluon connected Green function, represented by the full black blob in Fig. 1, and disconnected parts that are functions of the two gluon Green function, i.e. of the gluon propagator, the full blobs in grey.The four gluon connected Green function itself can be written in terms of 1-PI functions.In order to workout this decomposition let us consider the classical Yang-Mills field given by and introduce the generating functional for the one-particle irreducible correlation functions The one-particle irreducible Green functions are the functional derivatives of the Γ[A cl ] at vanishing classical fields.It follows from the above definitions that and, from this relation combined with Eqs ( 5) and ( 6), the orthogonality relation is derived.By taking functional derivatives of Eq (7), using the orthogonality relation (8), and after some straightforward algebra, one arrives at the equality whose diagrammatic representation is given in Fig. 2. The connected four gluon Green function is a sum of terms that include the four gluon 1-PI Green function and contributions proportional to the three gluon 1-PI Green functions.
On a lattice simulation it is the full Green function 4) and Fig. 1, that is measured.Its decomposition in terms of connected Green functions calls for the four-gluon connected Green function and disconnected terms that are proportional to products of gluon propagators.In the lattice formulation of QCD, the Green functions are proportional to the lattice volume and, therefore, the terms represented in Fig. 1 are such that the four gluon connected part is proportional to the lattice volume V , while the terms proportional to the propagators are proportional to V 2 .The lattice volume dependence of the terms that define G a1•••a4 µ1•••µ4 make it almost impossible to have a good Monte Carlo signal for the connected four gluon diagram, unless the contribution of the disconnected parts vanish.This can be achieved considering only those kinematical configurations such that the momenta p i and p j associated, respectively, with the external legs i and j are such that p i + p j ̸ = 0 for all i and j.From now on, it will be assumed that these conditions on the momenta of the external legs are fulfilled and, therefore, the contributions due to the disconnected diagrams can be discarded.
The four gluon connected Green function, represented diagrammatically in Fig. 2, is given by a term that includes the four gluon 1-PI Green function and diagrams that include three gluon 1-PI Green functions and gluon propagators.Then, to access the four gluon 1-PI Green function from the connected four gluon Green function, the three gluon contributions needs to be known.However, such information on the three gluon contributions is currently not available.One can either introduce some modelling for the three gluon 1-PI Green function, in such a way that its contribution can be removed, or further constrain the kinematical configuration of the external legs to suppress the contributions that call for three-gluon one-particle irreducible diagrams.For example, as discussed in [61], this is case when the incoming momenta are p 1 = p, p 2 = p, p 3 = p and p 4 = −3 p.A first and preliminary analysis of this kinematical configuration using lattice QCD simulations can be found in [66].
In general, for the Landau gauge, due to the orthogonality of the gluon propagator, for the kinematical configurations where the momenta of the external legs are all proportional, i.e. for the class of momenta such that where η and λ are real numbers, the terms in the decomposition of the four gluon connected Green function that require the contribution of the three-gluon vertex vanish.The proof of this statement can be found in App. A. Herein, in order to avoid any type of extra modelling, we choose to consider a subset of the kinematical configurations that belong to the class of momenta defined in Eq. (10).Although the modelling of the three gluon 1-PI is avoided, our choice based on the kinematics constraint the type of information that can be accessed in the simulation.
The choice of momenta is important to access the four gluon 1-PI Green function without extra assumptions.However, the color-Lorentz structure of this Green function is complex.For a general kinematical configuration, the number of tensors required to define a tensor basis to fully describe the four gluon 1-PI is large.For example, for the symmetric point, the four gluon vertex requires more than one hundred different tensors for its full description [63].General discussions of the tensor basis for this Green function can be found in [61,64].Ideally, one would access each of the form factors that multiply each tensor of the complete basis for the Green function with the help of a suitable projection operator.However, not only the lattice simulations convolute the 1-PI functions with the gluon propagators but also the quality of the signal-to-noise ratio in a Monte Carlo simulation favor certain types of projecting operators.In principle, the signal-to-noise ratio can always be improved by generating larger ensembles of gauge configurations, but a practical computation has always a limited statistical ensemble.
A complete tensor basis describing the four gluon 1-PI includes several types of operators that are proportional to the metric tensor and/or momenta.However, given that in a lattice simulation the measured function convolutes the four gluon 1-PI with gluon propagators, for the class of kinematical configurations considered here that is characterized by a single momentum, see Eq. ( 10), due to the orthogonality of the gluon propagator in the Landau gauge, the tensors proportional to momentum do not contribute to the full Green function that is measured in a simulation.This property simplifies considerably the tensor analysis the Green function and the tensors that contribute can only include the metric tensor for its Lorentz component and f abc or d abc , the fully anti-symmetric and fully symmetric structure constants, and/or δ ab in the description of the color part of the Green function.Then, by imposing Bose symmetry it is possible to build the possible tensor structures that, for the chosen kinematics, contribute to the Green measured on a simulation.A particular tensor that belongs to the class of allowed operators is the tensor structure that appears in the perturbative tree level four gluon Feynman rule Γ (0) abcd µνηζ = f abr f cdr (g µη g νζ − g µζ g νη ) + f acr f bdr (g µν g ηζ − g µζ g νη ) + f adr f bcr (g µν g ηζ − g µη g νζ ) .
From this tensor another operator that also contributes to the lattice Green function can be build replacing the fully anti-symmetric SU(3) structure constants f abc by the fully symmetric d abc and symmetrizing the Lorentz structure part.The corresponding tensor reads Γ (1) abcd µνηζ = d abr d cdr (g µη g νζ + g µζ g νη ) + d acr d bdr (g µζ g νη + g µν g ηζ ) + d adr d bcr (g µν g ηζ + g µη g νζ ) , Besides the two above tensors, one may also consider the tensor operator Γ (2) abcd µνηζ = δ ab δ cd + δ ac δ bd + δ ad δ bc g µν g ηζ + g µη g νζ + g µζ g νη (13) that was mentioned in [61].The three above tensors are not orthogonal to each other and their color-Lorentz scalars are where the numbers in the r.h.s.refer to the particular case where N = 3.In general, the above set of tensors do not define a complete tensor basis for the Green function that is measured on a simulation and others should be considered.However, in the current work only their contribution will be measured.
A. An orthogonal set of operators Although the operators Γ (0) , Γ (1) and Γ (2) are not orthogonal in the color-Lorentz space, by a linear combination of these operators it is possible to build three orthogonal operators from Γ (i) .Indeed, a straightforward calculation shows that the tensor operators t (3) abcd µνηζ = Γ (2) abcd µνηζ (22) are orthogonal in the color-Lorentz space and, therefore, where the normalization factors N j read once more the numbers in the r.h.s. are for N = 3.For a general kinematics, the three tensors t (1) , t (2) and t (3) do not define a complete basis neither for the full Green function nor for the one-particle function.The exception being the configuration of momenta such that p 1 = p 2 = p 3 = p and p 4 = −3 p where t (1) , t (2) and t (3) describe completely the correlation function [61].

B. Full Green function and lattice form factors
For the special kinematics under discussion herein, the full Green functions measured in lattice simulations gets only contributions from the connected four gluon full Green function.In the Landau gauge, the full Green function can be written as where F (p 2 1 , . . .), G(p 2 1 , . . .), H(p 2 1 , . . .), etc, are Lorentz scalar form factors, • • • represent the contribution of the remaining components of the tensor basis, supposedly but not necessarilly orthogonal to the space spanned by t (1) to t (3) , and is the orthogonal projector in momentum space that appears in the definition of the Landau gauge gluon propagator The 1-PI form factors to be measured from the full Green function are In general, the F (i) are linear combinations of the form factors F , G, H • • • that describe the 1-PI four gluon function.
It follows from the definitions ( 31) -( 33) that their r.h.s. are symmetric under interchange of any pair of momenta and Bose symmetry demands that the F (i) can only depend on the momenta through the combinations For the kinematical configurations investigated in the current work, where all momenta are proportional to each other and there is a unique momentum scale, one can write that F (i) ≡ F (i) (p2 ) to simplify the notation.
C. Lattice form factors F (i) and 1-PI form factors Our primer concern in this work is to compute the F (i) .These functions are not the form factors that describe the 1-PI four gluon vertex, see the definition in Eq. ( 28), but are given by linear combinations of F , G, H, . . .Assuming that the functions F , G and H give the dominant contributions to the form factors lattice form factors F (i) , then for the kinematical configurations (p 1 , p 2 , p 3 , p 4 ) = (0, p, p, − 2p) and (p 1 , p 2 , p 3 , p 4 ) = (0, p, 2p, − 3p), a straightforward calculation gives while for (p 1 , p 2 , p 3 , p 4 ) = (p, p, p, − 3p) one has If F and H are directly related to F (0) and F (2) , the form factor G(p 2 ) can be accessed by the linear combination of all the lattice form factors (41) We close this section by recalling the reader that for the kinematics (p 1 , p 2 , p 3 , p 4 ) = (p, p, p, −3p) the three tensor structures considered define a complete basis [61].

III. LATTICE SETUP AND CHOICE OF MOMENTA
The main goal of this work is to measure the form factors that describe the four gluon one-particle irreducible correlation function using lattice simulations.The results described below refer to simulations that use the Wilson action for β = 6.0.For this bare coupling constant the corresponding lattice spacing, measured from the string tension, is a = 0.1016(25) fm or, equivalently, 1/a = 1.943 (48) GeV.In the conversion from lattice to physical units we use the central values just quoted1 .The computer simulations were performed using the Chroma [68] and PFFT [69] libraries and run on the Navigator supercomputer at the University of Coimbra.
The form factors F (i) defined in Eqs (31) to (33) were computed with several ensembles of gauge configurations.Our first try was to consider the ensembles used in [5,70,71] that have 2000 gauge configurations for the 64 4 lattice and 1801 gauge configurations for the 80 4 lattice, both rotated to the Landau gauge 2 .For these ensembles of configurations the signal-to-noise ratios are quite poor, preventing from achieving any reliable measurement for the F (i) .To overcome the problem of the signal-to-noise ratio we have generated 9038 gauge configurations, in the Landau gauge, for a 32 4 lattice with the same β value, together with 4560 gauge configurations for a 48 4 lattice and β = 6.0.Our choice for β relies on the study of the gluon propagator [5] that suggests that for this β the lattice spacing effects are under control with possible mild volume effects in the IR that show up only in the smallest lattice.From now on, we will refer only to the measurements with these two latter lattices.
Our choice of momenta to measure the form factors F (i) considers only a subset of the class of momenta defined in (10).The subset of momenta was chosen with the aim to minimize the effects due to the breaking of rotational symmetry and our choice was to use momenta p i as close to each other as possible.The momentum configurations considered are In a simulation that uses hypercubic lattices, the breaking of rotational symmetry implies that the measured form factors F (i) are not only functions of p 2 but also of the H4 invariants p [n] [24, [72][73][74].If for the gluon propagator the building of these invariants and the evaluation of the effects due to the breaking of rotational symmetry is relatively straightforward, a similar analysis of Green functions with higher number of external legs becomes rather involved.
For the three point functions, the effects due to the breaking of rotational symmetry were discussed in [24] and for an hypercubic lattice, the breaking of rotational symmetry is minimized when one of the momenta is proportional to (1, 1, 1, 1).Similar conclusions have been observed for the two point functions [73][74][75] At least for the two point and three point gluon functions it was observed that the replacement of the naive momentum by an improved momentum, see Eq. ( 43) in Sec.IV, seem to handle the effects due to the breaking of rotation symmetry.
For the kinematical configurations under consideration, the form factors F (i) are functions of a single momentum scale.In order to improve the signal-to-noise ratio, given that the form factors are functions of p 2 , in all cases we average over equivalent momenta, including negative momenta, before performing the ensemble averages.For example, for each gauge configuration the form factors associated with the momenta are averaged before performing the ensemble averages.By doing so we are assuming that the dependence of the form factors on the H4 invariants p [n] is small or negligible.

IV. LATTICE FORM FACTORS
The continuum form factors F (i) defined in Eqs ( 31) to (33) were measured in the simulations using their lattice implementations as given in App.B. For the ensembles referred previously, it turns out that only the simulations using the 32 4 lattice have a good signal-to-noise ratio on a range of momenta.On the other hand, the results coming from the simulation with the 48 4 lattice show a good signal-to-noise ratio only for F (0) but not for the amputated functions.Anyway, despite the relative large errors associated with this latter lattice, the outcome of the simulation confirms the tendency observed in the 32 4 results and the two sets of data are compatible within one standard deviation.In all cases, the statistical errors were computed with the bootstrap method with a confidence level of 67.5%.The bare lattice data for the form factors will be described in terms of the improved lattice momenta where L is the length of the lattice.At least for the gluon propagator the use of the improved momenta handles part of the lattice effects.
The bare form factors F (0) , F (1) and F (2) associated with the full Green functions, i.e. that incorporate the contributions of the gluon propagators, are shown in Fig. 4 for all the kinematics and for all the momenta accessed in the simulation.The black circles are the data from simulations using the 32 4 lattice volume, while the red squares is the data measured on the 48 4 lattice volume.Note the different vertical scales that are associated with each of the FIG.3: Bare gluon propagator in the Landau gauge using the data from the simulations with the 32 4 and 48 4 lattices.The data reported is a subset of all the lattice momenta.Indeed, to handle the effects associated with the breaking of rotational symmetry we follow the recipe of [75] and perform the so-called cylindrical and conical cuts in momentum space, while for p < 0.7 GeV all data is included.
F (i) .Their order of magnitude can be understood looking at the bare gluon propagator function D(p 2 ), see Fig. 3, that is included in the data of Fig. 4.
The data reported shows there is a range of momenta where there is a clear sign for all the form factors that goes up to p ∼ 1 GeV.For higher momenta the statistical errors prevent the access to the behaviour of the various F (i) .Reduced errors can be accessed by using larger statistical ensembles of gauge configurations.However, the data also shows that the measurement of the four gluon Green functions is possible for gauge ensembles with, at least, about 10k configurations.
Let us now consider the same form factors but associated with the amputated Green functions.The form factors describing the 1-PI Green function, i.e. the amputated Green function, can be computed from the F (i) dividing the later functions by the gluon propagators that are associated with each of the external legs of the full Green function.For the size of each gauge ensemble, the statistical errors associated with the two point function are tiny, see the results of Fig. 3, and will be ignored when performing the division.Then, the statistical errors for the amputated Green function form factors come only from the statistical errors on the F (i) .
On a lattice simulation the (bare) computed functions relevant for our purpose are, using now a simplified notation, a 4 ⟨A(1) • • • A(4)⟩ for the four point function and a 2 ⟨A(1)A(2)⟩ for the propagator.From the division of the bare F (i) by the bare the gluon propagator, the would be amputated Green functions are To arrive at the corresponding dimensionless amputated Green functions, these ratios have to be multiplied by a 4 .Then, the bare lattice amputated Green function reads and can be seen in Fig. 5 as a function of the momenta.The data reported shows that there is a clear Monte Carlo signal for momenta up to ∼ 1 GeV for all the amputated form factors.For larger momenta the statistical errors become quite large and it is hard to comment on the behaviour of the amputated F (i) .The dependence on p shows an hierarchy between the various amputated form factors, with the form factor associated with tree level tensor structure F (0) being the largest of the three F (i) .In general, for the same form factor, the statistical error associated with kinematical configuration (p, p, p, −3p) is the largest.This can be understood looking at the lattice operators given in App.B. Indeed, for this kinematical configuration the number of averages performed for each configurations is In what concerns the dependence of the amputated form factors with momentum, the curves suggest that F (0) is essentially constant for p ≳ 0.6 GeV (although for p ≳ 1 GeV the large statistical errors prevent any firm conclusion) and decreases for smaller momenta.The systematics for F (1) and F (2) are harder to understand from Fig. 5.As discussed below, by taking into account the Bose symmetry and its implication on the amputated form factors gives a clear understanding on their dependence with momentum scales.
The gluon field is a bosonic field and, therefore, the form factors should depend on momenta as described in Eq. (34).In this sense, the form factors should be a function not of p but of and of the combination of scalar products of the various momenta mentioned in Eq. (34).Let us ignore this latter dependence and redo the plots, looking at the amputated form factors as a function of s and combining the various kinematics in the same plot.The lattice data for the three amputated form factors is reported in Fig. 6.In general, the lattice data for the amputated form factors is compatible with a dependence on s.For the level of statistical precision achieved in the simulation, the description of the lattice does not seem to require any further variable than s itself.The data of the amputated F (0) is compatible with a constant for s ≳ 1 GeV.However, F (0) is slightly suppressed at lower momenta.On the other hand F (1) and F (2) seem to be constant for s ≳ 1 GeV but increase when s becomes smaller.Moreover, the lattice data for the amputated F (1) suggests that this form factor should change sign at s ∼ 0.3 GeV.No change of sign is seen or suggested for F (2) .The relative strong increase in the F (1) and F (2)  observed when s approaches zero, seem to suggest that these form factors may have logarithmic divergences in the IR region.However, given the limited access to the low momenta region, the question of the logarithmic divergences is not solved within the current simulation.To answer this question one should have access to the deep infrared region that require simulations using larger lattices or with different lattice spacing.

V. SUMMARY AND CONCLUSIONS
Herein, the calculation of the four gluon one-particle irreducible Green function, in the Landau gauge, using lattice QCD simulations is addressed.The determination of the associated form factors from the full Green function is discussed and some of the 1-PI Green function form factors are computed.Our study shows that lattice simulations to access the four gluon 1-PI Green function require large ensembles of gauge configurations, meaning 10k or more gauge configurations per ensemble.Moreover, to resolve the 1-PI contribution to the full four gluon Green function, that is the Green function directly measured in the simulation, the kinematics of the external legs have to be carefully chosen.In particular, for the kinematics characterized by a unique momentum, it is possible to disentangle the four FIG.6: Dimensionless bare lattice amputated form factors F (i) as a function of s for the simulation on 32 4 lattice.The plots include the data for all the kinematical configurations.In the right plots the data for kinematical configuration with the large statistical error was omitted for a better reading on the dependence of the form factors with s. gluon 1-PI contribution from the disconnected parts and from the diagrams associated with the three gluon 1-PI in the full Green function.
In order to arrive at such large ensembles of configurations, the simulations reported use either a hypercubic lattice size of 32 4 , that has a physical volume of (∼ 3.3 fm) 4 , or a hypercubic lattice size of 48 4 , whose physical volume is about (∼ 4.9 fm) 4 .For the smaller lattice our ensemble has about 10K configurations, for the larger lattice the size of the ensemble is of the order of 4K.This prevents the larger lattice volume simulation to achieve a good signal-to-noise ratio and the statistical errors for the 1-PI form factors for the 48 4 lattice are quite large.Nevertheless, they confirm the tendency observed in the simulation with the 32 4 lattice.We are currently working on the production of larger statistical ensembles to improve our results, i.e. to achieve a good signal for larger range of momenta and to explore other kinematical configurations.The results of such simulations will be made public as soon as possible.FIG.7: Dimensionless bare lattice amputated form factors F (0) as a function of s for the simulation on 32 4 lattice as in Fig. 6 replacing the various data points for the same momentum with an weighted average using the statistical error squared as weight.
The complete description of the four gluon 1-PI irreducible Green function calls for a tensor basis with more than one hundred independent operators.However, for the kinematical configurations under consideration, the number of relevant operators that contribute to full Green functions is relatively small and, in some cases, can be handled exactly.Of the possible tensor operators that define the basis for the 1-PI, we have considered the three operators (11), ( 12) and ( 13), measured the associated form factors F (0) , F (1) , F (2) defined in ( 31), ( 32), (33), respectively, and their bare amputated versions that appear in the 1-PI four gluon irreducible Green function.
The tensor operator used to define the amputated form factor F (0) is the operator that appears in the tree level Feynman rule for the four gluon vertex.Our simulation gives an amputated form factor that is essentially constant for momentum above ∼ 0.5 GeV (see Fig. 6), with the data suggesting a suppression at small momenta.Despite the large ensemble, for momentum above ∼ 1 GeV it is difficult to disentangle its functional form.For momenta in the range ∼ 1 − 1.5 GeV, despite the size of the statistical errors, the data seems to be compatible with a constant value; this can be better seen by replacing the various data points for the same momentum by their weighted averages using the statistical error squared as a weight as reported in Fig. 7 for F (0) .This result is in qualitative agreement with the outcome of recent Dyson-Schwinger calculations [61,62,65] but not with the calculation described in [60].
The analysis of F (1) shows a negative form factor whose absolute value is about a quarter of the absolute value of F (0) , it has the opposite sign of this latter form factor and seems to approach zero as s is decreased.This behaviour suggests a change of sign for F (1) at s ∼ 0.15 GeV 2 and suggests a divergent behaviour for F (1) in the deep infrared region.Again, the observed F (1) is qualitatively in agreement with the Dyson-Schwinger result reported in [65].
The bare amputated F (2) form factor has, apart the possible change of sign, a similar behaviour as observed for F (1) , i.e.F (2) ∼ F (0) /4, with both form factors having the same sign and with the F (2) data suggesting a possible log divergence in the deep infrared region.The observed lattice data for F (2) is qualitatively in agreement with the recent Dyson-Schwinger calculation reported in [65].Note that the remaining Dyson-Schwinger studies do not access either F (1) or F (2) .The expressions in Eqs (B7) to (B15) will be used in the simulations to measure each of the lattice form factors F (i) .

FIG. 1 :
FIG.1: Diagramatic representation of the four gluon full Green function defined in Eq (4).The filled vertex in the left is the complete four point Green function, the blob filled in back is the connected Green function and the blobs in grey stand for one particle-irreducible diagrams.Color, Lorentz and momentum indices are omitted.

FIG. 2 :
FIG. 2: Diagramatic representation of the four gluon connected Green function (full black blob) appearing in Fig. 1 in terms of one particle-irreducible functions (full grey blobs).

FIG. 4 :
FIG. 4: Dimensionless bare lattice form factors F (i) for the various kinematics.Black dots are the data from the simulations on the 32 4 lattice, while the red squares refer to the simulation on 48 4 .Due to the large statistical errors the data from the larger volume is omitted in some of the plots.Note the different vertical scales used in each Fig. See main text for details.

FIG. 5 :
FIG.5: Dimensionless bare lattice amputated form factors F (i) for the simulation on 32 4 lattice.The data from the 484  simulation was omitted as it has a two large statistical error to add any useful information.For the lower momenta the data from the simulation with a 484 , that have smaller statistical errors, follow the trend observed in the 324 simulation and the two sets of data are compatible within one standard deviation.