Four-gluon vertex in collinear kinematics

To date, the four-gluon vertex is the least explored component of the QCD Lagrangian, mainly due to the vast proliferation of Lorentz and color structures required for its description. In this work we present a nonperturbative study of this vertex, based on the one-loop dressed Schwinger-Dyson equation obtained from the 4PI effective action. A vast simplification is brought about by resorting to ``collinear'' kinematics, where all momenta are parallel to each other, and by appealing to the charge conjugation symmetry in order to eliminate certain color structures. Out of the fifteen form factors that comprise the transversely-projected version of this vertex, two are singled out and studied in detail; the one associated with the classical tensorial structure is moderately suppressed in the infrared regime, while the other diverges logarithmically at the origin. Quite interestingly, both form factors display the property known as ``planar degeneracy'' at a rather high level of accuracy. With these results we construct an effective charge that quantifies the strength of the four-gluon interaction, and compare it with other vertex-derived charges from the gauge sector of QCD.

Instead, the nonperturbative aspects of the four-gluon vertex, to be denoted by IΓ abcd µνρσ , are relatively poorly known [23][24][25][76][77][78][79][80]; for perturbative studies, see [81][82][83][84][85][86][87][88].The main reason for this limitation is the large proliferation of Lorentz and color structures, which makes the treatment of this vertex exceedingly cumbersome in the continuum, and overly costly on the lattice.Since this vertex is inextricably connected with all other Green's functions through the coupled functional equations that govern their evolution, it is certainly desirable to improve our understanding of its dynamics.
Certain new insights gained from recent studies of the three-gluon vertex appear particularly promising for determining the leading nonperturbative features of the four-gluon vertex.In particular, the logistics of the three-gluon vertex are greatly simplified due to the property known as "planar degeneracy": with an error of less than 10%, the form factor associated with the classical tensor is the same for all momentum configurations lying on a given plane [58,89,90].Even though the dynamical reason that enforces this feature is not fully understood, its manifestation hinges crucially on the appropriate exploitation of the Bose symmetry of the three-gluon vertex, a symmetry shared also by the four-gluon vertex.It is, therefore, possible that the planar degeneracy may be a property of the fourgluon vertex, leading to a great simplification of practical computations.Parallel to these developments, lattice simulations have been carrying out exploratory studies of the four-gluon vertex in simplified kinematics [91][92][93], and are expected to access a wider array of momentum configurations in the near future.
In the present work we study the quenched four-gluon vertex (no dynamical quarks) in the context of the one-loop SDE derived from the four-particle irreducible (4PI) effective action [94][95][96][97][98] at four loops [99,100].Note that, within this formalism, the extremization of the n-loop n-PI effective action with respect to the m-point correlation functions (m ≤ n) generates the corresponding equations of motions (SDEs), which are expressed in terms of (n − m + 1)-loop diagrams [96,99,100].We emphasize that our analysis does not treat the entire system of coupled SDEs that govern the propagators and vertices entering in the SDE of the four-gluon vertex: instead, we consider the four-gluon SDE in isolation, using lattice results as inputs for all other dynamical components.
The focal point of our attention is the transversely projected four-gluon vertex, IΓ abcd µνρσ , which is precisely the one simulated on the lattice in the Landau gauge [91][92][93], and appears in the majority of physical applications [80,101,102].The analysis is restricted to the case of the "collinear" kinematics, where all vertex momenta are parallel to each other, or, equivalently, proportional to a single momentum p (i.e., p i = x i p).This kinematic choice, together with the proper exploitation of the charge conjugation symmetry, drastically reduces the tensorial structure of IΓ abcd µνρσ : only 15 elements, t abcd i,µνρσ , are required for the full description of its Lorentz and color content, with the associated form factors denoted by G i (i = 1, . . ., 15).Note that this basis contains as one of its elements the transversely projected classical (tree-level) tensor of the four-gluon vertex; in particular, t 1 = Γ 0 .
There is an additional simplification that we will incorporate in our study, prompted by the following two observations.First, the elements t 1 , t 2 , and t 3 satisfy a special orthogonality condition with respect to the rest of the basis.Second, they depend only on p, being independent of the variables x i .These facts reduce significantly the algebraic complexity of the problem, facilitating the determination of the attendant form factors G 1 , G 2 , and G 3 ; we therefore focus exclusively on this special subset of contributions.
In the 4PI formalism, all propagators and vertices appearing in the one-loop diagrams comprising the SDE that governs the evolution of the G i are fully dressed, including the four-gluon vertices.This converts the SDE into an integral equation, involving the unknown function IΓ abcd µνρσ evaluated in different kinematic configurations.The approximation used for dealing with this complication is to assume that the tree-level form factor G 1 is dominant, and that it satisfies planar degeneracy, in exact analogy to its three-gluon vertex counterpart.
The results of our study may be summarized as follows: (i) Under the approximations mentioned above, the form factor G 2 vanishes identically.
(ii) G 1 and G 3 are computed for several collinear configurations.We find that G 1 is suppressed in the infrared with respect to its tree-level value, G (0) 1 = 1, reaching a finite value at the origin, while G 3 displays a logarithmic divergence in the deep infrared, originating from the ghost loops.
(iii) For momenta below 1 GeV, our results for G 1 are compatible with planar degeneracy to within 1.7%.In the case of G 3 , a similar tendency is observed, albeit with a lesser degree of accuracy.
(iv) With the help of the G 1 computed for an extensive set of collinear configurations, we construct a "band" of effective charges, α 4g (p), which may be used to quantify the strength of the four-gluon interaction.The α 4g (p) is then compared with the effective charges α 3g (p) and α cg (p), obtained from the soft-gluon limits of the three-and ghost-gluon vertices, respectively.For momenta below 1.5 GeV, we find the clear hierarchy α 3g (p) < α 4g (p) < α cg (p), in agreement with earlier results presented in [8,25].
The article is organized as follows.In Sec.II we summarize the general features of IΓ abcd µνρσ in collinear configurations.In Sec.III we show how the charge conjugation symmetry prohibits the presence of color tensors of the type f abx d cdx .In Sec.IV we define the Lorentz tensor basis for collinear configurations, and determine the projectors that allow us to extract the attendant form factors.In Sec.V we present the SDE governing the evolution of IΓ abcd µνρσ , derived from the 4PI effective action at the four-loop level, and its renormalization.In Sec.VI, the SDE is solved for general collinear configurations, and the form factors G 1,3 are determined and analyzed.In Section VII we compare our results to previous continuum studies, and comment on the impact of the three-gluon vertex in the four-gluon SDE.In Sec.VIII we present an analytic proof of the ultraviolet finiteness of the form factors G 1,3 , which complements our numerical findings.Next, in Sec.IX, we compute the four-gluon effective charge.Finally, our conclusions are discussed in Sec.X, and certain technical points are summarized in two Appendices.

II. FOUR-GLUON VERTEX IN COLLINEAR KINEMATICS
The starting point is the definition of the four-gluon correlation function as the vacuum expectation value of the time-ordered product of four SU(3) gauge fields, A a µ (p), in momentum space, with p 1 +p 2 +p 3 +p 4 = 0.The diagrammatic representation of G abcd µνρσ , shown in the upper panel of Fig. 1, distinguishes between connected and disconnected contributions.The connected Green's function, denoted by C abcd µνρσ , has the general form where C abcd µνρσ is the amputated vertex.C abcd µνρσ may be further separated into one-particle irreducible (1PI) and one-particle reducible (1PR) contributions, namely where ∆ ab µν (q) = −iδ ab ∆ µν (q) denotes the gluon propagator, and IΓ abc αβγ (q, r, p) the full threegluon vertex [see lower panel of Fig. 1].
In general kinematics, the vertices IΓ abcd µνρσ and IΓ abcd µνρσ have a proliferation of tensorial and color structures, which give rise to a large number of form factors.In the present work we restrict ourselves to the nonperturbative study of IΓ abcd µνρσ in collinear configurations.This type of configurations are defined through a single four-momentum p, and all vertex momenta p i (i = 1, 2, 3, 4) are proportional to it, namely where the parameters x i are real, and satisfy x 1 + x 2 + x 3 + x 4 = 0 due to four-momentum conservation.
The reason for this particular kinematic choice is twofold.First, the tensorial structure of the four-gluon vertex simplifies enormously, rendering the resulting dynamical equations completely tractable.Second, it allows for the direct comparison with contemporary lattice simulations [91][92][93], where these configurations are employed as benchmarks in the initial stages.
Note, in fact, that, while the SDEs access directly the vertex IΓ abcd µνρσ or IΓ abcd µνρσ , the lattice computes the full G abcd µνρσ .When all momenta are collinear, the disconnected contributions may be eliminated simply by imposing the additional restriction x i + x j ̸ = 0, which makes propagator-like transitions impossible, due to the induced momentum non-conservation.
Alternatively, appropriate projectors may be chosen [see Sec.IV], which remove the disconnected contributions from the simulation without restricting the values of the x i .As for the 1PI terms, when Eq. (2.10) is fulfilled, all three-gluon vertices in the 1PR diagrams are of the form IΓ αβγ (p) = A(p 2 )(p α g βγ + p β g αγ + p γ g αβ ) + B(p 2 )p α p β p γ , and therefore IΓ αβγ (p) = 0. Thus, the term in the curly bracket of Eq. (2.7) vanishes, and one finally isolates IΓ abcd µνρσ on the lattice [78].

III. COLOR STRUCTURE AND CHARGE CONJUGATION SYMMETRY
It is well-known that the color structure of the four-gluon vertex is comprised by the appropriate combinations of three basic elements, namely f abc , the totally symmetric d abc , and the identity matrix δ ab .The 15 combinations built out of these elements are [78,81] and their permutations.The number of possibilities is reduced to 9 thanks to 6 identities, valid for general SU(N ) [78,81], namely and two permutations for each.Moreover, for N = 3, the additional identity [81] reduces the number of independent color tensors to 82 .
It turns out that the invariance of the theory under charge conjugation prohibits the presence of terms of the type f d, thus eliminating 3 additional color combinations.In earlier works [8,76,103], and under special assumptions, such terms have been shown to vanish; however, no exact proof justifying their complete omission has been presented in the literature.
The general implications of charge conjugation symmetry for the Green's functions of a pure SU(N ) theory were worked out in [104].Specifically, a Green's function with n color indices, Γ a 1 a 2 ...an (we suppress Lorentz indices and momenta), must satisfy In the above equation, A ab is a certain color matrix that conjugates the charges; its concrete form will not be necessary here.It suffices to state the following properties of A ab : {1} A ab ̸ ∝ δ ab , i.e., A ab is not proportional to the diagonal matrix.
{2} A is an orthogonal matrix, i.e., A ⊺ = A −1 , where ⊺ denotes matrix transposition; To see how Eq. (3.4) can then be used to constrain the color structures of vertex functions, consider first the simpler case of the three-gluon vertex, IΓ abc αµν (q, r, p).The most general structure possible is given by [81] IΓ abc αµν (q, r, p) = f abc V αµν (q, r, p) + d abc U αµν (q, r, p) , where the tensors V αµν (q, r, p) and U αµν (q, r, p) may be decomposed in an appropriate basis.
Using the Jacobi identity and the second line of Eq. (3.2), it is easy to verify that the expression in Eq. (3.5) satisfies the first line of Eq. (3.4).Then, substituting Eq. (3.5) into the second line of Eq. (3.4) and using properties {3} and {4}, we obtain Hence, charge conjugation symmetry implies that the three-gluon vertex is proportional to Let us now consider the four-gluon vertex.From property {2} follows that and thus Hence, multiplying by A aa ′ A bb ′ A cc ′ A dd ′ , we find where we used property {3} to obtain the last equality.
Similarly, using property {2}, we get Finally, using property {4}, it is straightforward to show that Consequently, the number of independent color structures in the four-gluon vertex is reduced to 5, namely f abx f cdx ; f acx f bdx ; δ ab δ cd ; δ ac δ bd ; δ ad δ bc . (3.12)

IV. TENSOR BASIS AND PROJECTORS
In this section we discuss the tensorial basis that will be used for describing the four-gluon vertex in collinear configurations, and the projectors that allow us to extract particular form factors.
From this point on, we will employ the short-hand notation to indicate a quantity A (with Lorentz and color indices, as well as scalar form factors) in general collinear kinematics.In the case of a specific configuration, i.e., (1, 1, 1, −3), we will be reverting to the explicit notation, i.e., A(p, p, p, −3p).
First, notice that in collinear configuration one can form 10 Lorentz tensors with four indices and one independent momentum [25,85], covering all linearly independent combinations of the forms g µν g ρσ (3 tensors) ; g µν p ρ p σ (6 tensors) ; Combining the 5 independent color tensors of Eq. (3.12) with the Lorentz tensors of Eq. ( 4.3), we see that for collinear configurations IΓ abcd µνρσ is comprised by 15 color-Lorentz tensors.However, the basis resulting from the tensor product of the building blocks in Eqs.(3.12) and (4.3) is not manifestly Bose symmetric.
To amend this, we use linear combinations of the above elements to construct a basis of 15 new tensors t abcd i,µνρσ , and associated form factors Notice that in the decomposition given by Eq. ( 4.4), one has that (i) Each t i is manifestly Bose symmetric.Consequently, the form factors G i are individually symmetric under the exchange of any two components x i implicit in x.
(ii) The t i , and therefore the G i , are all dimensionless.This property prevents the form factors G i from developing kinematic divergences when one of the momenta vanishes.
Moreover, the different form factors can be compared directly.
(iii) Finally, the tree-level tensor Γ 0 defined in Eq. (2.12) is an element of the basis, namely To construct a basis satisfying the above properties, we employ the S 4 permutation group formalism developed in [105].Specifically, we combine the multiplets of momenta, colors, and Lorentz indices computed there into S 4 singlets (invariants).The details are presented in Appendix A, and the expressions for the resulting t i are given in Eq. (A13).
In order to extract the form factors G k from the SDE, we construct projectors P k that isolate specific elements of the basis.In compact notation, we have where the symbol "⊙" denotes the full contraction of all tensor and color indices, namely To determine the projectors P k , we first contract both sides of Eq. (4.4) with an arbitrary t j , where we define the symmetric matrix C with entries [C] ji = (t j ⊙ t i ).Next we multiply both sides of Eq. (4.6) by the inverse with k = 1, . . ., 15.
The first three tensors, t 1,2,3 , of the basis of Eq. (4.4) form a rather special subset.To begin, they are orthogonal to the remaining tensors, i.e., In fact, t 3 , is orthogonal to all other t i .Furthermore, we find that t 1,2,3 are the only tensors of the basis that are independent of x i3 .Indeed, for t 1,2,3 the expressions in Eq. (A13) reduce to the relatively compact forms These two observations have important consequences for the projectors P 1 , P 2 and P 3 which extract the associated form factors.The orthogonality expressed in Eq. (4.8) gives rise to a significant simplification in the procedure of Eq. (4.7), allowing these projectors to be expressed as linear combinations of only t 1,2,3 ; in particular, P 3 is simply proportional to t 3 .As a consequence, we find that P 1 , P 2 , and P 3 , are themselves independent of x 1 , x 2 and x 3 ; through the procedure outlined, we find that they read We hasten to emphasize that, since we made explicit use of Eq. (3.3), the above projectors are valid for SU(3), and may not hold for a general SU(N ) group.
In what follows we will concentrate on the form factors G 1,2,3 (x, p), since they do not mix with the remaining form factors in any collinear configuration.
We conclude by pointing out that the projector P 1 completely eliminates disconnected contributions when applied to G abcd µνρσ , because its contraction with color structures of the type δ ab δ cd vanish.Thus, in principle, the tree-level tensor structure can be isolated on the lattice, even for x i + x j = 0, by employing P 1 .On the other hand, the tensors t 2,3 mix with disconnected terms in G abcd µνρσ , since the projectors P 2,3 do not annihilate δ ab δ cd .

V. SDE OF THE FOUR-GLUON VERTEX
In this section, we set up the SDE that governs the form factors G 1,2,3 of the fourgluon vertex for an arbitrary collinear configuration, and discuss the renormalization of the resulting equations.
We employ the version of the four-gluon SDE derived from the formalism of the 4PI effective action [94][95][96][97][98] at the four-loop level [99,100]; we remind the reader that, within this formalism, the SDE of a given Green's function is obtained by extremizing the variation of the effective action with respect to this particular function.
It turns out that the SDE of the four-gluon vertex derived from the 4PI effective action is automatically symmetric with respect to permutations of its external legs, and contains fully dressed propagators and vertices as its ingredients, depicted diagrammatically in Fig. 2. The four-gluon SDE built out of these components is shown in Fig. 3; diagram (d 1 ) is accompanied by five additional permutations, (d 2 ) by two, (d 3 ) by five, and (d 4 ) by two [78].In what follows we will denote by the sum of each representative diagram and its permutations.
To obtain the SDE of the transversely-projected four-gluon vertex, IΓ abcd µνρσ (p 1 , p 2 , p 3 , p 4 ), one simply contracts both sides of the SDE in Fig. 3 by the transverse projectors corresponding to the four external legs; the additional projectors, needed for the internal lines, + permutations come automatically from the gluon propagators, which in the Landau gauge, assumes the completely transverse form, ∆ µν (q) = ∆(q 2 )P µν (q).Thus, in addition to IΓ abcd µνρσ , the various diagrams depend on the transversely-projected three-gluon and ghost-gluon vertices,
Thus, the SDE for IΓ abcd µνρσ may be expressed in the following compact way where the ( ds i ) denotes the transversely-projected counterpart of the (d s i ) defined in Eq. (5.1).We now turn to the renormalization of the SDE given in Eq. (5.3).All quantities appearing in that equation are bare (unrenormalized); the transition to the renormalized quantities is implemented multiplicatively, by means of the standard relations where the subscript "R" denotes renormalized quantities, and Z A , Z c , Z 1 , Z 3 , Z 4 , and Z g are the corresponding (cutoff-dependent) renormalization constants.In addition, we employ the exact relations which are imposed by the fundamental Slavnov-Taylor identities (STIs) [106,107].Substituting the relations of Eq. (5.7) into Eq.( 5.3) and using Eq.(5.8), it is straightforward to derive the renormalized version of Eq. (5.6), given by where the subscript "R" in ds i R abcd µνρσ denotes that the expressions given in Eq. (5.3) have been substituted by their renormalized counterparts.We emphasize that the Z 4 multiplied by the tree-level structure is the only renormalization constant to be determined.Note in particular that, because all vertices are fully dressed, none of the renormalization constants appears multiplying any of the contributions ds i R abcd µνρσ ; thus, the usual complications associated with multiplicative renormalization are absent [99,100].
In order to determine Z 4 we adopt a variant of the momentum subtraction (MOM) scheme [108,109].Within this scheme, the gluon and ghost propagators assume their treelevel values at the subtraction point µ, i.e., In general studies of vertices, these two conditions are supplemented by an additional one, fixing the value of a form factor (or combination thereof) at some special kinematic configuration [110][111][112].For example, in the recent SDE analysis of [90], the renormalization of the three-gluon vertex was such that the leading form factor would acquire its tree-level value in the soft-gluon configuration [see Eq. (3.13) therein].
For the case of the four-gluon vertex we employ an analogous scheme.In particular, we single out the configuration (0, p, p, −2p), and require that There is nothing special about this configuration, other than its relative numerical simplicity; any other choice would be equally good for renormalizing the SDE.
Note that, when the scheme defined by Eqs.(5.10) and (5.11) is employed, the conditions used for renormalizing the three-gluon and ghost-gluon vertices may not be simultaneously employed, because this would lead to a violation of the relations of Eq. (5.8) imposed by the STIs [113].The renormalized three-gluon and ghost-gluon vertices used in our calculations as inputs of the four-gluon SDE must undergo finite renormalizations, which will adjust their values such that the relations of Eq. (5.8) are fulfilled.This procedure of finite renormalization is described in detail in Appendix B.
The determination of Z 4 proceeds by first employing the special projection of Eq. (4.10) to define where x 0 denotes the reference configuration (0, 1, 1, −2) in Eq. (5.11).Then, imposing the condition Eq. (5.11) at the level of Eq. (5.9), we find Note that all contributions to d 0 1 (p 2 ) originating from the ghost boxes vanish, i.e., P 1 ⊙ ( ds 1R )(x, p) Then, substituting Eq. (5.13) into Eq.(5.9), we obtain the renormalized SDE for IΓ When no ambiguity can arise we drop the index "R" to avoid notational clutter.
The final renormalized expressions for both G 1,3 (x, p) are obtained by acting on Eq. (5.14) with the associated projectors given in Eq. (4.10), The SDEs in Eq. ( 5.15) will be solved in the next section, under a set of well-founded assumptions, which will be discussed there in detail.It turns out that, under these assumptions, one may demonstrate analytically the ultraviolet finiteness of the resulting G i (x, p); the detailed proof is presented in Sec.VIII.

VI. RESULTS: EMERGENCE OF PLANAR DEGENERACY
In this section we solve the SDE of the four-gluon vertex under certain simplifying assumptions, and check the extent of validity of the planar degeneracy at the level of the solutions obtained.

A. Preliminary considerations
For the SDE treatment, we employ appropriate inputs for the various functions entering in the kernels of Eqs.(6.7) and (6.8).In particular, for the gluon propagator, ∆(r 2 ), and the ghost dressing function, F (r 2 ), we use fits given by Eqs.(C11) and (C6) of [114], respectively, to the lattice results of [42,54].
It is clear that the arguments of the G i appearing in the integrands of the SDE are comprised by combinations x i p and the loop momentum k, a fact that complicates the numerical treatment.However, a considerable simplification occurs if (i ) the form factor G 1 is assumed to be dominant, making negligible all dependence on other form-factors, in particular G 2,3 , and (ii ) the G 1 is assumed to display "planar degeneracy" at some notable level of accuracy.In such a case, one implements the substitution and where, in general, The form of G * 1 (s 2 ) will be determined by means of an iterative procedure, outlined in the next subsection.Note that, for a given G * 1 (s 2 ), all four-gluon vertices in the integrands of Eq. (5.15) will be replaced by Eqs.(6.4) and (6.5), and evaluated at the value of s2 corresponding to each vertex.For example, the s2 assigned to the vertex IΓ admn µσαβ (p 1 , p 4 , k, −k 3 ) in the diagram d3 in Eq. ( 5.3) is given by s2 = 1 2 (p 2 1 + p 2 4 + k 2 + k 2 3 ), where k 3 := k+p 1 +p 4 .There are two main consequences stemming from the implementation of these two approximations.
First, the form factor G 2 (x, p) vanishes identically.Indeed, we find that the diagrams ( d1,2 ) of Eq. ( 5.3) cannot generate the particular color combination associated with G 2 [cf.Eq. (4.9)], independently of the approximation used for the ghost-gluon and three-gluon vertices.On the other hand, the diagrams ( d3,4 ) can contribute to G 2 , but only if non-treelevel tensor structures of the four-gluon vertices therein are dressed.Since a nonzero G 2 can only arise from the effect of dressing non-tree level tensor structures, or including higher loop diagrams in the SDE, it is reasonable to expect it to be subleading in comparison to G 1,3 .In particular, it is evident that, perturbatively, G 2 can only be nonzero at two loops or higher.
Second, the remaining equations for G 1 (x, p) and G 3 (x, p) decouple, assuming the schematic form where the K i are the various kernels defined from the diagrams ds i (x, p); all functions other than G * 1 have been absorbed into them.The term quadratic in G * 1 (s 2 ) stems from graph (d 4 ); and x 0 denotes the configuration (0, 1, 1, −2) in Eq. (5.11).
In order to proceed with the numerical solution, the Eqs.(6.7) and (6.8) must be passed to Euclidean space, following standard conventions (see, e.g., Eq. (5.1) of [119]).For the numerical treatment, the resulting integrals are expressed in hyper-spherical coordinates.
Numerical integration is performed with the Gauss-Kronrod implementation of [120].
We integrate over a grid of external squared momenta distributed logarithmically over the interval [10 −4 , 10 4 ] GeV 2 .configurations, defined in (6.9).G 1 (x, p) was computed employing the approximation IΓ = Γ 0 on the r.h.s. of the Eq.(6.7).The blue curve on the right panel represents G av 1 (s 2 ), the average of these four configurations in terms of s2 , defined in Eq. (6.11).

B. Planar degeneracy as an infrared feature.
The first step in determining G * 1 (s 2 ) is to set G * 1 (s 2 ) = 1 everywhere on the r.h.s. of Eq. (6.7).After this substitution, any collinear configurations for G 1 (x, p) are obtained from Eq. (6.7) by simply carrying out the corresponding integrations.Thus, we determine numerically G 1 (x, p) in four collinear kinematic configurations, whose coordinates (x 1 , x 2 , x 3 , x 4 ) are given by (1, 1, 1, −3) , (0, 1, 1, −2) , (0, 1, 2, −3) , (1, 2, 3, −6) .(6.9) This choice is random; there is nothing particular about these four cases.The results are displayed in Fig. 4: on the left panel, where the four configurations are plotted as functions of the momentum p, they lead to clearly distinct curves; however, when plotted in terms of the Bose symmetric variable s of Eq. (6.6), they merge into each other in the range 0 ≤ s ⪅ 1 GeV, manifesting clearly the onset of the planar degeneracy in this region of momenta.Note that in the case of collinear momenta, Eq. (6.6) reduces to This promising result motivates the iteration of the procedure, in order to obtain the optimal form of G * 1 (s 2 ).In particular, we determine the numerical average, G av 1 (s 2 ), of the four configurations defined in Eq. (6.9), i.e., G 1 (x j , p) ; (6.11) the result of this "averaging" is shown as the blue continuous curve on the right panel of Fig. 4.Then, the second step in the iterative procedure is taken, by setting on the r.h.s. of Eq. (6.7) G * 1 (s 2 ) = G av 1 (s 2 ) ; (6.12) in this way, the G av 1 (s 2 ) acts now as the new "seed" for obtaining the next generation of results for the same four configurations.These new results will then be averaged and fed back into the r.h.s. of the original set of four equations.This iterative procedure concludes when the relative error between two consecutive averaged solutions is smaller than 0.1%; typically, convergence is achieved after about 15 iterations.In Fig. 5, we compare the G * 1 (s 2 ) corresponding to the first (blue dashed curve), and last iterations (red continuous curve).
The final result is clearly suppressed with respect to the initial G * 1 (s 2 ) = 1 (tree-level).Note that this particular procedure decouples effectively Eq. (6.7) from Eq. (6.8): the form factor G 3 (x, p) may be obtained by means of a single integration, once the final G * 1 (s 2 ) has been computed from Eq. (6.7) and substituted into Eq.(6.8).
However, even though G 3 (x, p) is not subject to an iterative procedure, it is advantageous to define a curve G * 3 (s 2 ), which potentially approximates the general collinear G 3 (x, p), and allows this form factor to be reconstructed from a single curve (see Fig. 7).Moreover, the curve G * 3 (s 2 ) will serve as a reference for quantifying how well G 3 (x, p) satisfies the property of planar degeneracy.Therefore, by analogy with Eqs.(6.11) and (6.12), we set where the sum is over the same kinematic configurations used to define G * 1 (s 2 ), given in Eq. (6.9).The resulting G * 3 (s 2 ) is shown on the right panel of Fig. 5 for the case G * 1 (s 2 ) = 1 (blue dot-dashed) and for the final G * 1 (s 2 ) (red continuous).As is clear from Fig. 5, for s ≥ 0.3 GeV, G * 3 (s 2 ) is considerably smaller than G * 1 (s 2 ).On the other hand, for s → 0, the form factor G * 3 (s 2 ) diverges logarithmically, whereas G * 1 (0) is finite.As shown in [78], the divergence of G * 3 (s 2 ) results from the massless ghost loop (d 1 ) of Fig. 3, which contributes to G 3 (x, p), but not to G 1 (x, p).The rate of divergence can be determined by an asymptotic analysis along the lines of the one performed for the three-gluon vertex in Appendix A of [90].Specifically, we obtain where Z 1 is the renormalization constant of the ghost-gluon vertex.Note that Z 1 is finite in the Landau gauge [106,109,121]; within our renormalization scheme, it turns out to be very close to unity, Z 1 = 1.01 .
In order to test the reliability of the entire procedure, we compute G 1 and G 3 for fourteen random collinear kinematics, using as input in Eqs.(6.7) and (6.8) the G * 1 (s 2 ) obtained at the end of the iteration procedure (red continuous curve in Fig. 5).Ten of these configurations are obtained by employing a random number generator for the x i , while the last four correspond to those used earlier for the determination of G * 1 (s 2 ); the coordinates (x 1 , x 2 , x 3 , x 4 ) of all of them are listed below4 : In Fig. 6 we plot both form factors G 1,3 (x, p), for the configurations listed in Eq. (6.15), as functions of the momentum p (left panels) and the Bose symmetric variable s (right panels).
The remarkable agreement between the fourteen curves when they are plotted as functions of s is to be contrasted to the considerable disparity observed when they are plotted as functions of p. Thus, for s < 1 GeV, the approximation appears to be particularly robust.
To quantify the level of accuracy of Eq. (6.16) within the range 0 < s < 2 GeV, we compute the percentage error defined as where i = 1, 3, and G i (s 2 j ) denotes any one of the fourteen curves shown on the right panels of Fig. 6, while G * i (s 2 ) are our reference curves, shown as red continuous lines in Fig. 5. Starting with the classical form factor, G 1 (x, p), we find that for s ≤ 1 GeV, the error δ 1 (s 2 j ) ≤ 1.7%.As can be seen already from Fig. 6, the relative separation between curves of different kinematics grows as s increases; consequently, in the range 1 ≤ s ≤ 2 GeV, δ 1 (s 2 j ) increases, reaching the maximum value of 4.4%.
For G 3 (x, p), notice that, for all the kinematics shown in Fig. 5, the curves display zero crossings at around s = 0.7 GeV, where the relative error δ 3 (s 2 j ) becomes ill-defined.Away from the crossing, we find that δ 3 (s 2 j ) ≤ 7.0% for s ≤ 0.4 GeV, whereas for 1 ≤ s ≤ 2 GeV the δ 3 (s 2 j ) increases, reaching the value δ 3 (s 2 j ) ≤ 33.7%; however, in this latter range of momenta, G * 3 (s 2 ) is itself rather small.In addition to the inspection of discrete sets of kinematic configurations, one may visualize the accuracy of the planar degeneracy approximation more generally through 3D plots, where the x i components vary continuously.
In Fig. 7, we plot G 1,3 (x, p) (yellow surfaces) as functions of the momenta x 1 p and x 2 p.
To make this visualization possible, we reduce the number of independent variables by fixing p 3 = (x 1 + x 2 )p.In the same figure, we map G * 1,3 (s 2 ) into the (x 1 p, x 2 p) plane, by setting p 3 = (x 1 + x 2 )p in Eq. (6.10), such that s2 = 3(x 1 p) 2 + 3(x 2 p) 2 + 5(x 1 p)(x 2 p).With this mapping, the curves for G * i (s 2 ) (red continuous in Fig. 5) become surfaces, shown in blue in Fig. 7.For both form factors, G 1,3 (x, p), the blue and yellow surfaces show a high degree of overlap in the infrared, in agreement with Eq. (6.16), while differences become visible with increasing momenta.
Let us finally point out that, if planar degeneracy were an exact property of the G i (x, p), there should be no dependence on x for fixed s; instead, for an approximate planar degeneracy, the G i (x, p) should show a mild dependence on x.The residual dependences of the G i (x, p) on x can be visualized by plotting surfaces corresponding to fixed values of s; their deviation from the absolute flatness will be then a measure of the exactness of the planar degeneracy.
To this end, first note that for a given vector x, the value of p that gives rise to a given value of s can be uniquely determined by inverting Eq. (6.10).Specifically,  where x 4 has been eliminated by using momentum conservation.Hence, we can determine how G 1 (x, p) varies as a function of x 1 , x 2 and x 3 , for fixed s: from an available array of values for G 1 (x, p), we pick the one value that corresponds to the unique p obtained from Eq. (6.18), once an s and a set of (x 1 , x 2 , x 3 ) have been substituted in it.
In Fig. 8 we first set x 1 = 1, and then plot G 1,3 (x, p) for fixed values of s and general values of x 2 and x 3 .Specifically, we show surfaces corresponding to s = 0.5 GeV (brown), s = 1 GeV (blue), and s = 2 GeV (green).In both panels, it is clear that for s = 0.5 GeV and s = 1 GeV, the G 1,3 (x, p) are almost perfectly constant, i.e., nearly independent of x 2 and x 3 .Hence, we confirm once more the validity of the planar degeneracy property expressed by Eqs.(6.5) and (6.16) in the infrared.
For s = 2 GeV, the surfaces of G 1,3 (x, p) in Fig. 8 begin to deviate from perfect flatness, in line with our previous observations that planar degeneracy breaks down with increasing s.However, even in this case, the surfaces of G 1,3 (x, p) are rather flat for most values of x 2 and x 3 .The only exceptions occur near the borders of the plot, when at least one of the x i is small, with particularly stronger dependences on x when both x 2 , x 3 ≈ 0.

VII. COMPARISON WITH PREVIOUS WORKS: IMPACT OF THE THREE-GLUON VERTEX
In this section, we compare our results with previous continuum studies of the four-gluon vertex.We recall that in the present work we employed a 4PI-based truncation of the SDE, given by Eq. ( 3), where all vertices in the one-loop diagrams appear dressed.On the other hand, in [78] all vertices on the r.h.s. of the SDE are set to tree-level, whereas in [8,25] one vertex per diagram appears at tree-level.As such, the comparison of these results to ours allows us to assess the effect that the dressing of the vertices, and especially of the three-gluon vertex, has on the four-gluon SDE.
For this task, we will focus on the configuration (p, p, p, −3p), corresponding to c 11 of Fig. 6, which has also been studied in [8,25,78]; for simplicity, we limit our discussion to the classical form factor G 1 .
Note that [25,78] employed renormalization schemes different to the one used in the present work, defined by the prescription in Eqs.(5.10) and (5.11).Hence, to perform a meaningful comparison, we rescale all results to be equal to 1 at µ = 4.3 GeV, i.e., In Fig. 9 we show our result for G 1 (p, p, p, −3p) as a continuous red line, while the result of [78] is shown in blue dot-dashed, and that of [25] corresponds to the green dotted curve.
The most prominent difference is seen for p < 1.5 GeV.In this region, the results of [25,78] display a large peak, which is nearly absent in our G 1 (p, p, p, −3p).Indeed, in this range of momenta, our result has a local maximum of 0.9 at p = 0.8 GeV.In contrast, the result of [78] is strongly enhanced in the infrared, reaching a maximum of 1.5 at 0.4 GeV, while [25] obtains a smaller maximum of 1.1 at 0.6 GeV, and becomes very similar to ours for p ⪅ 0.4 GeV.
To understand the origin of these differences in the infrared behavior, we have computed the individual contributions of each diagram of Fig. 3, and compared to the Fig. 5 of [78] and Fig. 8 of [25].We find that, in the infrared, the diagrams (d 2 ) and (d 3 ) contribute positively to G 1 (p, p, p, −3p), in agreement with [25,78], but vary considerably in size within the different truncations.
Then, we note that both of these diagrams contain three-gluon vertices.As has been firmly established by lattice [43, 49-51, 112, 116, 117] and continuum studies [58,98,115,[122][123][124], the three-gluon vertex is strongly suppressed in the infrared in comparison to its tree-level value.Hence, in our treatment, the diagrams (d 2 ) and (d 3 ) yield minimal contributions, due to the suppressing effect of dressing the three-gluon form factor L sg (s 2 ) [cf.Eq. ( 6.1)].Consequently, the shape of our G 1 is dominated by the tree-level and the negative (d 4 ) diagram.
In contrast, in [78], where all vertices on the r.h.s. of the SDE are at tree-level, (d 2 ) and (d 3 ) give large contributions, giving rise to the peak seen in Fig. 9. Finally, the truncation of [25] yields a result between ours and that of [78], since there is only one tree-level vertex per diagram.
From these observations we conclude that the approximation employed for the three-gluon vertex has a profound impact on SDE analyses of the four-gluon vertex.
Finally, we point out that the four-gluon vertex has also been studied in the context of scaling solutions in [8,25,77].In this case, the classical form factor displays infrared divergences that are absent in our framework.

VIII. ULTRAVIOLET FINITENESS OF THE FORM FACTORS G 1 AND G 3
Even though the numerical evaluation of the SDEs presented in the previous section gave rise to results that are free of ultraviolet divergences, it would be clearly advantageous to devise an analytic demonstration of this fact.In this section, we show that when the SDE undergoes the renormalization procedure outlined in Sec.V, and under the assumptions discussed in Sec.VI A, one obtains ultraviolet finite form factors G 1 and G 3 .In particular, G 1 becomes finite after a single subtraction, while G 3 turns out to be automatically free of ultraviolet divergences, exactly as expected.
We begin by writing the SDEs in Eq. (5.15) in Euclidean space using hyperspherical coordinates; then, integrating over the angle between k and p, and setting y := k 2 , we arrive at where Λ 2 UV is an ultraviolet cutoff, and f i (y, p 2 , x) denotes the sum of the integrands of the projected diagrams.The closed expressions for the f i (x, p) are rather complicated; however, for the purposes of renormalization, only their behavior for asymptotic y is relevant.
At this point we introduce an arbitrary mass scale M , with p 2 , µ 2 ≪ M 2 , and split the integrals in Eq. ( 8.1) into a finite and a potentially divergent part, i.e., The type of divergence contained in the second integral may be established by performing a Taylor expansion of f i (y, p 2 , x) around p 2 /M 2 = 0, where the ellipsis denotes higher derivatives.Note that, upon setting p 2 = 0, all reference to the parameter x := (x 1 , x 2 , x 3 , x 4 ) drops out from f i (y).
Evidently, for the renormalization to go through, the only potential divergence must be isolated in the term f i (y), whereas the terms with the derivatives must yield finite integrals when substituted into Eq.(8.2).
To check if this indeed so, let us first focus on f 3 (y); a detailed calculation reveals that it consists of four contributions, I i (y), one from each diagram ( ds i ) [i = 1, 2, 3, 4], namely where I 1 (y) = y −1 B 4 1 (y, y, π)F 4 (y) ; I 2 (y) = y 3 ∆ 4 (y)L 4 sg (y) ; To further evaluate the I i (y), we assume that, in the limit y → ∞, the functions comprising them reproduce their one-loop resummed forms [8,19,121,125] i.e., where c = 1/ ln(µ 2 /Λ 2 MOM ), and Λ MOM is the QCD mass-scale in the MOM scheme.Note that this assumption is built into the fits employed for all quantities that serve as inputs for the SDE: indeed, the expressions for F (y), ∆(y), and L sg (y) given respectively in Eqs.(C6), (C11) and (C12) of [114] are designed to capture the correct one-loop resummed anomalous dimension, while the numerical results in Fig. 13 of [58] satisfy B 1 (k, −k, 0) → 1 at large k.
The anomalous dimension of G 1 can be obtained by standard renormalization group analysis [126] from the one-loop form of the renormalization constant Z 4 , given in Eq. (B6).
Turning to the form factor G 1 , note that only ( ds 3 ) and ( ds 4 ) contribute to it, yielding which upon integration yields the divergent contribution of G 1 , given by the constant Clearly, the subtraction prescribed by the first relation in Eq. (5.15) will cancel the divergence in Eq. (8.12), giving rise to a finite G 1 .Specifically, the exact same steps leading from Eq. ( 8.1) to Eq. (8.12) may be repeated for f 1 (y, µ 2 , x 0 ), expanding around µ 2 /M 2 = 0; x 0 denotes the reference configuration (0, 1, 1, −2) in Eq. (5.11).
We end this section by confirming that the higher derivative terms in Eq. ( 8.3) are ultraviolet finite, and do not interfere with renormalization.In particular, with the aid of Eq. (8.6), we obtain for constants a i , b i and c i that depend in general on x, and the ellipsis denotes terms that are subleading for large y.Integrating the r.h.s. of Eq. (8.13) over y yields a convergent result; keeping higher order derivatives in the expansion of Eq. ( 8.3) leads to even faster converging integrals.
It is important to emphasize that, for asymptotic values of the momentum, our numerical solutions for G * 1 are completely compatible with the anomalous dimension of 2/11, quoted in Eq. (8.6); this is shown in Fig 10.

IX. FOUR-GLUON EFFECTIVE CHARGE
Typically, the quantitative description of the strength of a given interaction proceeds through the construction of the corresponding renormalization group invariant (RGI) ef-fective charge [70,77,[127][128][129][130][131][132][133][134].In the case of the four-gluon interaction, the appropriate RGI combination is obtained with the aid of the relation Z Specifically, introducing the gluon dressing function Z(p 2 ) := p 2 ∆(p 2 ), it is straightforward to establish that the combination retains the same form before and after renormalization, and thus, does not depend on the renormalization point µ, event though the individual components are µ-dependent.Then, the quantity α 4g (x, p) will be constructed using inputs renormalized at a given µ, which we choose µ = 4.3 GeV; evidently, any other choice of µ would be equally good, since the value of α R s (µ 2 ) is adjusted to compensate precisely the µ-dependence in the other components.In order to obtain a concrete realization of a vertex-derived effective charge, one normally specifies a particular kinematic configuration for the vertex form factor entering in it.For instance, in the case of the effective charges α cg (p) and α 3g (p), obtained from the ghost-gluon and three-gluon vertices, respectively [70,77,[128][129][130], where B sg (p 2 ) := B 1 (p, −p, 0) and L sg (p 2 ) are the "soft-gluon" limits of the classical form factors of the ghost-gluon and three-gluon vertices, respectively.Instead, in order to obtain a more representative picture for the strength of the four-gluon interaction, we evaluate α 4g (x, p) using not one, but all fourteen collinear configurations listed in Eq. (6.15).Thus, we generate a sequence of effective charges, which end up comprising the blue band shown in Fig. 11.
For comparison, in the same figure we show also the ghost-gluon and three-gluon effective charges of Eq. (9.2).For momenta p ⪅ 1.5 GeV there is a clear hierarchy, which is expressed by the inequality α 3g (p) < α 4g (x, p) < α cg (p), independently of the collinear configuration chosen for α 4g (x, p).In fact, the same hierarchy was obtained in [8,25] for charges defined by the three-and four-gluon vertices in their totally symmetric configurations, and the soft-ghost kinematics for the ghost-gluon vertex.
In the ultraviolet, the different charges slowly merge into each order.Nevertheless, even at p = 5 GeV, the α 4g , α 3g and α cg are visibly different.Moreover, for p ⪆ 3.5 GeV, the charge hierarchy is modified with respect to the infrared: α 4g (x, p) < α cg (p) < α 3g (p).To verify that this ordering and the overall sizes of the charges in the ultraviolet are correct, rather than a numerical or truncation artifact, we compute the relation between these charges at one loop.
Perturbatively, one charge can be expanded in powers of any other [113]; in particular, where the ellipses denote higher powers of α 4g .Then, our task reduces to determining the one-loop coefficients a 1 and b 1 .Using the results of [135] for B sg (p 2 ), F (p 2 ), Z(p 2 ), and , together with the one-loop result for G 1 (x 0 , p) given in Eq. (B7), we find that Since a 1 > b 1 > 0, Eq. ( 9.3) implies that the hierarchy found in Fig. 11 at large p agrees with perturbation theory.In fact, it is straightforward to confirm that the agreement is not just qualitative, but quantitative.In particular, setting p = 4.3 GeV and α 4g = 0.23 [see Eq. (B11)] into Eq.( 9.3), we get α cg = 0.24, and α 3g = 0.26; these values match those obtained from the full nonperturbative calculation (i.e., the corresponding curves in Fig. 11), namely α cg = 0.26 and α 3g = 0.27, within 8% and 3%, respectively.

X. CONCLUSIONS
We have presented a comprehensive nonperturbative study of the transversely-projected four-gluon vertex by means of the SDE obtained within the 4PI effective action formalism.
For the purposes of this work we have restricted our attention to the case of the collinear kinematics, which simplify considerably the complexity associated with this vertex.This kinematic choice, in conjunction with the elimination of color terms of the type f abx d cdx by virtue of the charge conjugation symmetry, reduces the number of possible form factors down to fifteen.A special subset of three form factors, including the one associated with the tree-level (classical) tensorial structure, is subsequently singled out and studied in detail.
Our treatment of the SDE uses as external inputs results obtained from lattice QCD and previous continuous studies.
The most outstanding nonperturbative feature that emerges from our analysis is the "planar degeneracy" of the form factors considered.Specifically, at a high level of accuracy, all kinematic configurations that lie on a given plane share the same form factors in the infrared.
This special property, first observed at the level of the three-gluon vertex [58,89,90,115], appears to be shared by the four-gluon vertex, at least in the context of the collinear configurations that we have studied.It would be clearly important to probe the extent of validity of the planar degeneracy for more complicated configurations of momenta, such as general soft-gluon, where one momentum vanishes and the other three are completely arbitrary; calculations in this direction are already in progress.
It would be of course particularly interesting to confront our findings with lattice QCD, as has been the case with all other fundamental vertices of the theory.Unfortunately, the preliminary data available [91][92][93] are not sufficiently conclusive for confirming any of the main features encountered in the present study.It would be interesting to refine the lattice simulations, in order for a meaningful comparison with the results of functional methods to become feasible.In this context, the assumption of planar degeneracy might prove particularly useful, because it would increase significantly the available statistics.
Let us finally point out that in our analysis we have not considered dynamical quarks.
The inclusion of quarks may proceed straightforwardly, by making two basic modifications to the treatment presented here.First, in the diagrammatic representation of the SDE in Fig. 3, one has to include additional diagrams containing quark-loops [136,137].Second, denote generic doublets and triplets, respectively.Then we define the dot ("•") product as D • D ′ = aa ′ + ss ′ , T • T ′ = uu ′ + vv ′ + ww ′ , (A9) together with "star" ("⋆") and "wedge" ("∧") products given by are S 4 singlets [105].Hence, to construct a Bose symmetric tensor basis for collinear configurations, we form all possible combinations of the above types using the multiplets in Eqs.(A3) through (A7), and select from these 15 linearly independent singlets.
From the procedure outlined above, we obtain the following basis C ) .
Note that the dependence on the x i appears only through the scalar multiplets, s, D and T , of Eqs.(A3) and (A4).Hence, the tensors t 1,2,3 are independent of x i , whereas all others depend on it.
Note that for the gluon and ghost propagators the renormalization prescriptions of Eqs.(B1) and (5.10) coincide.Therefore, the multiplicative renormalizability expressed by Eq. (5.7) implies that the renormalization constants Z A and Z c have the same values in both schemes, i.e., where the index "B" denotes unrenormalized quantities 6 .It follows that ∆ asym (q 2 ) = ∆(q 2 ) and F asym (q 2 ) = F (q 2 ).
On the other hand, the vertex form factors and the value of the coupling are different in the two schemes considered.Specifically, Eq. (5.7) implies L sg (r 2 ) = Z 3 Z asym

3
L asym sg (r 2 ) , B 1 (q, r, p) = Then, it is easily shown from the STIs of Eq. (5.8) that Z 3 = (Z 4 Z A ) 1/2 , Z 1 = Z 3 Z c Z −1 A , and , and identical relations with "asym" superscripts in all quantities.Hence, using Eq.(B2), the ratios in Eq. (B3) can be expressed as In what follows, we assume that the renormalization point µ = 4.3 GeV is sufficiently large such that Z 4 /Z asym 4 can be approximated by a one-loop calculation.
We begin by computing the one-loop Z 4 in the renormalization scheme of Eq. (5.11).To that end, we substitute in Eq. (6.7) the tree-level forms of all inputs, i.e., ∆(q 2 ) → 1/q 2 ,

. 11 )
Combining the transformation properties given by Eqs.(3.9), (3.10) and (3.11) into the second line of Eq. (3.4), one can easily show that the form factors of the color structures of the form f d must all vanish.

Figure 3 .
Figure 3. Diagrammatic representation of the SDE for the full four-gluon vertex, IΓ abcd µνρσ , derived from the 4PI effective action at the four-loop level.Contributions obtained through the permutations of the external momenta in the various diagrams are omitted.

Figure 4 .
Figure 4. Comparison of G 1 (x, p) expressed as a function of the momentum p (left panel) and as a function of the Bose symmetric variable s (right panel), for the four collinear kinematic

Figure 6 .
Figure 6.The form factors G 1,3 (x, p) for the fourteen random collinear kinematic configurations, defined in Eq. (6.15), plotted as functions of momentum p (left panels) and as functions of the Bose symmetric variable s (right panels).The disordered pattern seen on the left panels is to be contrasted with the tight clustering of the curves displayed on the right ones.

Figure 7 .
Figure 7. Form factors G 1,3 (x, p) (yellow) compared to the approximation G * 1,3 (s 2 ) (blue), mapped into the (x 1 p, x 2 p) plane using Eq.(6.10), and setting x 3 = x 1 + x 2 .In both cases, the advantage of planar degeneracy becomes evident: entire surfaces are reconstructed, to a great accuracy, from the knowledge of a single curve.

Figure 10 .
Figure 10.Comparison of the nonperturbative solution for G * 1 (y) with its one-loop resummed behavior given by Eq. (8.6).

4 /Z 4 .
order to convert between the asymmetric MOM and the prescription of Eqs.(5.10) and (5.11), we simply need to obtain the finite, and regulator independent, ratio Z asym