A Dyson-Schwinger study of the four-gluon vertex

We present a self-consistent calculation of the four-gluon vertex of Landau gauge Yang--Mills theory from a truncated Dyson--Schwinger equation. The equation contains the leading diagrams in the ultraviolet and is solved using as the only input results for lower Green functions from previous Dyson--Schwinger calculations that are in good agreement with lattice data. All quantities are therefore fixed and no higher Green functions enter within this truncation. Our self-consistent solution resolves the full momentum dependence of the vertex but is limited to the tree-level tensor structure at the moment. Calculations of selected dressing functions for other tensor structures from this solution are used to exemplify that they are suppressed compared to the tree-level structure except for possible logarithmic enhancements in the deep infrared. Our results furthermore allow one to extract a qualitative fit for the vertex and a running coupling.


Introduction
The non-perturbative analysis of quantum field theories is one of the great challenges in physics. One particular approach is to use functional equations for Green functions which are the basic building blocks of a quantum field theory. In quantum chromodynamics (QCD) they can be used as input for hadron phenomenology and strong-interaction matter studies, see, e.g., [1][2][3][4][5][6][7][8][9][10][11], but they also allow direct conclusions on non-perturbative aspects like confinement [12][13][14] or dynamical mass generation, see, e.g., [1,2,6] and the references therein.
The main challenge for functional methods is the necessity to truncate the originally infinite hierarchy of functional equations for Green functions and to quantify the resulting uncertainties. The most straightforward way to assess the quality of a particular truncation by going beyond it in a systematic way as provided, e.g., by derivative or vertex expansions is often rather difficult. Hence alternative possibilities for tests are welcome and widely used, such as comparisons with results from lattice simulations where they are available.
In this paper we focus on Yang-Mills theory in the Landau gauge. The good understanding of this particular covariant gauge that was established in the past provides the basis for many of the more phenomenological investigations of QCD. The Landau gauge propagators have been well studied with various methods, e.g., lattice simulations [15][16][17][18][19][20], Dyson-Schwinger equations (DSEs) [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37], the functional renormalization group (FRG) [32,38], a variational approach [39], a one-loop model calculation with gluon mass term [40], or the (refined) Gribov-Zwanziger framework [41][42][43][44][45][46]. Also three-point functions are by now better understood [47][48][49][50][51][52][53][54][55][56][57][58] and their equations of motion can be solved self-consistently [53,57,58]. The qualitative behavior of propagators and vertices is well-captured by standard truncations of functional equations, but their quantitative reliability still needs to be tested and improved. Based on our most recent results for the complete set of two and three-point functions [57], however, there is quite compelling evidence to expect that the system of DSEs truncated to the primitively divergent Green functions might yield a rather good approximation that does not rely on any further external input. The two pieces missing to confirm this are the four-gluon vertex, which was the only remaining model input in such calculations [57,58], and the two-loop diagrams in the gluon propagator DSE (of which one also contains the four-gluon vertex). Few direct calculations of the latter exist [59][60][61], but for the four-gluon vertex information is even more scarce. Here we provide further information from calculating the four-gluon vertex within a state-of-the-art truncation that takes into account its full momentum dependence. An interesting additional feature of this truncation for the four-gluon vertex DSE is that for the first time no model input is required here, since only primitively divergent lower n-point Green functions enter at this level which are all known sufficiently well.
So far, little non-perturbative information on the four-gluon vertex is available to compare to, unfortunately. Even for the three-gluon vertex, available lattice results are limited to very restricted kinematical regions [48]. Due to the existence of six kinematic variables in the four-gluon vertex (as compared to three for three-point functions) the situation is even much more difficult here. Thus, even if lattice data for the four-gluon vertex will become available in the future, a kinematically reasonably complete coverage will likely remain impossible for some time to come. A continuum method has a clear advantage in this respect, although the resolution of the full kinematic dependence is certainly a challenge here as well. Perturbative results at the symmetric point were presented in Ref. [62,63]. Studies beyond perturbation theory can be found in Refs. [64][65][66]. In Ref. [65] the box diagrams were studied in a certain momentum configuration which we will refer to as configuration A below. As input non-perturbative propagators from the so-called scaling type were used [21,22,27]. In Ref. [66] this was extended to include all UV leading one-loop diagrams and propagators of the decoupling type were used as input [67].
In this work we go beyond these previous studies in several ways. First of all, we take into account the full momentum dependence. This is useful when our results are used as input in future calculations, because they provide a guideline to develop approximations that still capture the main features but are easier to handle than the full results. The full momentum dependence is also required to solve the equation self-consistently so that we can study the back coupling effects for the vertex. Finally, we use for the first time non-perturbative input for the three-gluon vertex that is in good agreement with lattice data.
In Sec. 2 we fix our notations and present a self-contained derivation of the four-gluon vertex DSE. This section also contains information on the truncation, the tensor basis, the kinematics and the renormalization. The input we employ is described in Sec. 3 and our results are presented in Sec. 4. We conclude in Sec. 5. Some technical details about color contractions, tensor bases and the numerical calculations to solve the four-gluon vertex DSE can be found in the appendices.

Derivation of the four-gluon vertex DSE
The Lagrangian density of Yang-Mills theory, fixed to the linear covariant gauge, is

2)
3) where A is the gluon field and (c) c is the (anti-)ghost field. The gauge fixing parameter is denoted by ξ. For the Landau gauge it will be set to 0 later. The field strength tensor F a µν and the covariant derivative D ab µ are given by From the Lagrangian density the path integral is defined as where J, σ andσ are the sources for the gluon and ghost fields. For the derivation of Dyson-Schwinger equations the one-particle irreducible (1PI) action will be used that is obtained from the path integral via a Legendre transformation: Here J i represents the sources J, σ andσ and Φ ∈ {A cl , c cl ,c cl } denotes the classical fields determined by In the following we will drop the subscript cl again. For convenience we use a multi index for the fields and sources that includes field species, Lorentz and color indices and position (or alternatively momenta). Consequently, repeated indices include summation over the discrete and integration over the continuous variables. For example, ). The derivation of DSEs starts from the integral over the derivative of the path integral whichassuming no boundary terms exist -must vanish, viz.
where S = dxL eff YM is the gauge-fixed action. With Eq. (2.8) this is rewritten to the master equation for the DSEs of 1PI Green functions, see, e.g., [1,68,69]: This means that one should differentiate the action S with respect to a field φ i and then replace every field φ i by the classical field Φ i plus a functional derivative multiplied by D J ji . The DSEs resulting from this procedure can conveniently be represented by Feynman diagrams. D J ji is the second derivative given by If we set the sources to zero, this becomes the propagator D J=0 ji . In our case we have the ghost and the gluon propagators, where P T is the transverse projector, P T µν (p) = g µν − p µ p ν /p 2 . From the master equation (2.11) the DSE for any n-point function can be derived by n − 1 field derivatives and subsequently setting the sources to zero. For n = 2 one obtains the DSEs of the inverse propagators and for n > 2 the equations for the vertices which we define as 1 As a characteristic feature every diagram in a DSE contains a bare vertex, which entails restrictions on the set of possible diagrams. Due to the self-interaction of the gluons via the three-and the four-gluon vertices, the DSE of the four-gluon vertex is the four-point function with most terms, namely 60 terms of which there are 20 one-loop and 39 two-loop diagrams. In particular the tensorial structure of the four-gluon vertex is, because of its four color and four Lorentz indices, considerably more complicated than those of the pure ghost or the ghost-gluon four-point functions. So in principle the derivation of its DSE is straightforward but tedious, and we used the Mathematica package DoFun [70,71] for this task. Figure 1. The truncated four-gluon vertex DSE. We use the following shorthand notation: {i, j, k} represents three diagrams where the indices at the first, second and third positions are chosen. We call the second and the third diagram on the right-hand side swordfish (SF) and dynamic triangle (DT), respectively. The diagrams in the second line are named gluon box (GlB), ghost box (GhB) and gluon triangle (GlT) (from left to right). Further, we call the swordfish and the dynamic triangle the dynamic diagrams since they depend on the four-gluon vertex. Accordingly, we call the diagrams of the second line static diagrams.
With the truncation discussed in Sec. 2.2, the DSE for the four-gluon vertex Γ abcd µνρσ (p, q, r, s) is schematically written as Γ abcd µνρσ (p, q, r, s) = Γ (0),abcd µνρσ + Λ abcd µνρσ (p, q, r, s) + Λ acdb µρσν (p, r, s, q) + Λ adbc µσνρ (p, s, q, r) + . . . , (2.16) where we used that all one-loop diagrams appear in three permuted versions. Suppressing momentum arguments, the sum of all unpermuted one-loop diagrams, denoted by Λ abcd µνρσ in Eq. (2.16), is (see also Fig. 1) The subscripts denote the names of the diagrams as explained in Fig. 1. The five individual i Λ abcd µνρσ are called primitive diagrams from which other diagrams are constructed by permutation. They are given by Purely gluonic vertices are denoted by Γ with the corresponding number of Lorentz and color indices, while the ghost-gluon vertex has only one Lorentz index which is associated to the first color index.
Here the renormalization constants Z 1 , Z 1 and Z 4 of the three-gluon, the ghost-gluon and the fourgluon vertices appear.

Asymptotic behavior and truncation
Our truncation scheme consists of two parts: As usual, we discard several diagrams based on their asymptotic behavior. As it happens, all remaining Green functions required for the calculation are already known and we need no model input. However, we further simplify the system by a restriction of the color and Lorentz bases for the four-gluon vertex. The latter aspect of the truncation is discussed in Sec. 2.3.
For the diagrammatic truncation of the four-gluon vertex DSE we follow the same arguments as employed for the three-gluon vertex [57,58]. The main guidelines are the correct IR behavior and the inclusion of all diagrams contributing at one-loop order in the UV. Thus we retain only one-loop diagrams with primitively divergent Green functions. To be precise, from the set of all 20 one-loop diagrams we discard one with a ghost-gluon five-point function, one with a gluonic five-point function and three ghost triangles with a ghost-gluon four-point function leaving the 15 diagrams representing the Λ's in Eqs. (2.16) and (2.17) which are depicted in Fig. 1.
For the discussion of the IR behavior of the vertex we have to elaborate shortly on the IR behavior of Landau gauge Yang-Mills theory in general. It is known that the system of propagators allows two different types of solutions [29,31,32]. One is called decoupling solution and actually consists of a family of solutions that have all in common that the ghost dressing function G(p 2 ) and the gluon propagator Z(p 2 )/p 2 become constant in the IR: This solution is, besides in the DSE approach [29][30][31][32]36], also found by many other methods like lattice calculations, e.g., [15][16][17][18][19][20], with the functional renormalization group [32], within the refined Gribov-Zwanziger framework [44][45][46], with an effective model [40] or in a variational approach [39]. From the functional perspective the different solutions are distinguished by the boundary condition imposed for the ghost propagator DSE [32]. As a limiting case also the solution c → ∞ exists. This is called the scaling solution and characterized by a power law behavior of all Green functions [21,33,72]. The propagator dressing functions can be written as The exponents δ gh and δ gl are related by 2δ gh + δ gl = 0. Typically they are given in terms of κ := −δ gh . Its value can be calculated analytically as κ = (93 − √ 1201)/98 ≈ 0.6 [24,25]. This value can only change if the ghost-gluon vertex is not regular in the IR, viz. if its IR limit depends on the angle between two momenta [25,73]. Within modern truncations for the ghost-gluon vertex no such dependence was seen [53].
A very convenient feature of the scaling solution is that the qualitative behavior of all Green functions can be derived without truncations by combining the two systems of functional equations given by the FRG and the DSEs [33]. For a vertex with 2n ghost and m gluon legs the dressing functions behave qualitatively as (p 2 ) (n−m)κ where p is an IR momentum scale [33,72]. Consequently we expect that for this type of solution the four-gluon vertex behaves as (p 2 ) −4κ . This behavior is exhibited by all diagrams with a bare ghost-gluon vertex. Within our truncation these are the ghost boxes. However, other diagrams with ghost-ghost-gluon-gluon or ghost-ghost-gluon-gluon-gluon functions exist that have by power counting the same IR behavior but are discarded here. In fact, all diagrams can be classified in terms of their scaling behavior as determined by the bare vertex.
For the decoupling type of solution such a straightforward classification is not possible. From the three-gluon vertex it is known that it diverges logarithmically in the IR [54,[56][57][58] and it was conjectured that this is also true for the four-gluon vertex [56]. For both quantities these divergences stem from the ghost loops. We will come back to this point in Sec. 4.

Tensor basis
The four-gluon vertex is undoubtedly the most complicated primitively divergent Green function of Landau gauge Yang-Mills theory. With four color indices it possesses a non-trivial color structure and the four Lorentz indices allow a multitude of Lorentz tensors. We start with a discussion of the former.
As specified in Sec. 3, we use only the totally anti-symmetric structure constant f abc for the threepoint functions. To our knowledge no proof exists that the totally symmetric color part is non-zero in three-point functions. Thus we neglect the totally symmetric d abc from the beginning, but note that the color tensors we use can partly be also expressed via the d symbols. The building blocks are then the Kronecker delta δ ab and the totally anti-symmetric structure constant f abc . The basis constructed from them is (2.21) Another possible tensor, C abcd

6
= f adn f bcn , is not included since it can be expressed via the Jacobi identity as C abcd . Furthermore, contractions of more anti-symmetric structure constants can be reduced to this set. In particular, in the four-gluon vertex DSE terms of the form and for SU (2) to For SU (N ) with N > 3 the tensor C abcd 7 is linearly independent and must be considered as well. This is shown in Appendix A.
The Lorentz space is even more complicated than the color space. If one constructs all possible Lorentz tensors from the metric tensor δ µν and the three independent momenta, one arrives at 138 tensors. They can be split into the following classes: 3 dimensionless tensors: δ µν δ ρσ , δ µρ δ νσ and δ µσ δ νρ 54 tensors of dim. 2: However, from considerations along the lines of Refs. [58,74,75], it turns out that there are only 136 independent tensors [76]. In Landau gauge the completely transverse subspace is sufficient, which still contains 43 linearly independent tensors [64].
For a first study of the vertex within our truncation the full transverse basis is still by far too large. Thus we restrict ourselves here to the tree-level structure of the four-gluon vertex which is given by and replace all full four-gluon vertices by The non-perturbative information is contained in the dressing function D 4g (p, q, r, s). This strategy worked very well for the three-gluon vertex. We will comment in Sec. 4 on its applicability for the four-gluon vertex.

Bose symmetry
A truncated DSE in general does no longer reflect all the symmetries of the full equation. For gluonic Green functions the Bose symmetry is of special importance. If we simply went ahead and calculated the truncated DSE Eq. (2.16), the results would not be symmetric under the exchange of the leg attached to the bare vertices and another one. This would entail that the results depend on the way the four-gluon vertex is fed back into the DSE. Thus it is necessary to symmetrize the results. We want to stress that this effect comes from using dressed vertices in our setup and would be absent if all vertices were bare. The straightforward way for symmetrization is to average over the four DSEs with different momenta at the legs attached to the bare vertices. Within our truncation this corresponds to calculating all possible permutations of the primitive diagrams. The actual number of diagrams can be reduced using inherent symmetries of the equation. Nevertheless, this would lead to an increase in complexity and computing time, which we will avoid as explained below. First of all, we consider how to extract the dressing function D 4g (p, q, r, s) from the DSE given in Eq. (2.16). For this the following projection is employed that explicitly gets rid of all longitudinal parts (momentum arguments are suppressed on the right-hand side): Note that Λ adbc µσνρ does not represent the sum of all one-loop terms of the four-gluon vertex DSE but appears in three permuted versions in the DSE, see Eq. (2.16). From Eq. (2.27) one can obtain the symmetrized, transversely projected tree-level dressing function by D 4g (p, q, r, s) = Z 4 + 1 4 L(p, q, r, s) + L(q, r, s, p) + L(r, s, p, q) + L(s, p, q, r) L(p, r, s, q) + L(q, s, p, r) + L(r, p, q, s) + L(s, q, r, p) L(p, s, q, r) + L(q, p, r, s) + L(r, q, s, p) + L(s, r, p, q) . (2.28) Visualization s r q q r s s r q Table 1. Definitions of the kinematic configurations A, B and C used in plots. The dashed lines represent the fourth momentum vector given by momentum conservation.
Equation (2.28) is obtained from Eq. (2.16) by projecting it onto the transverse tree-level structure and symmetrizing it. Naively, one would expect 4! = 24 terms. However, this number reduces to 12 since some diagrams turn out to be identical. The reasons are the Bose symmetry of the four-gluon vertex itself and the irrelevance of the direction of the loop momentum, see Fig. 1.
To reduce the computational effort we did not calculate all L explicitly. Rather, we computed the normalized one-loop expression, given by Eq. (2.27), and from that the dressing function with Eq. (2.28). In other words, calculating L(p, q, r, s) with full momentum dependence gives us access to all other variants of L appearing in Eq. (2.28) so that D 4g (p, q, r, s) can be computed from the calculation of only one L.
We explicitly tested what happens when no symmetrization is employed and found that there is a considerable impact on the results. Furthermore, we want to mention that no transverse projection was employed in Ref. [65]. Thus our results for the ghost and gluon boxes cannot directly be compared to theirs. Unfortunately, the transverse projection also increases the complexity of the four-gluon DSE significantly, by about an order of magnitude.

Kinematics
Another new aspect of our investigation is that we take into account the full momentum dependence of the vertex. This is due to the large number of kinematic invariants a considerably complex task. The vertex depends on three independent momenta, say, s, r and q, from which six kinematic invariants can be formed, e.g., s 2 , r 2 , q 2 , s · r, s · q and r · q. However, this choice has the disadvantage that the domains of the latter three are not independent. Hence, we directly use spherical coordinates to describe the three momenta at the cost of having some sets of coordinates that describe the same momentum vectors. 2 Exploiting the O(4) symmetry we define where S, R, Q ∈ R + and θ r , θ q , ψ q ∈ [0, π] . From Eq. (2.29) it is easy to see that the domain of q · r is not [−QR, QR] but depends on the other angles. For example, for θ r = π/2 we have q ·r = QR sin(θ q ) cos(ψ q ). This can in general not be rewritten into a form QR cos(α r ) with α r ∈ [0, π].
As a second example consider θ r = 0. Then q · r = QR cos(θ q ), but θ q is already the free angle from s · q. Thus we prefer to work with the angle ψ q instead of r · q.
The dressing function itself is defined on a six-dimensional grid of the variables s 2 , r 2 , q 2 , θ r , θ q and ψ q . Typically we use 15 points for the radial and 7 points for the angular variables. When we have a converged solution also points on a finer grid are calculated.
For visualization of our results we have to fix some variables. In Tab. 1 three kinematic configurations are shown that will be used later. Configuration A is the one also used in Refs. [65,66].

Renormalization
The integrals of the four-gluon vertex DSE are logarithmically divergent. However, since we are using input that was obtained within the MiniMOM scheme [21,77], we are not free to subtract these divergences via a momentum subtraction. Within that scheme the renormalization constant of the ghost-gluon vertex Z 1 is fixed to 1 for the Landau gauge. The ghost and gluon propagators, on the other hand, were obtained from a self-consistent calculation which also entails certain values for their renormalization constants Z 3 and Z 3 , respectively. The Slavnov-Taylor identities (STIs) then fix the renormalization constants of the three-and four-gluon vertices as The corresponding values can be found in Table 2.
We use the given value for Z 4 in the tree-level expression of the vertex. However, the renormalization constants in front of the integrals in Eq. (2.18) are treated differently. Motivated by the gluon propagator equation, where this is necessary to obtain the correct anomalous dimension, we replace the renormalization constants Z 1 and Z 4 by momentum dependent functions [22,53]. The purpose of these functions is to effectively add a UV dressing to the bare vertices which brings the equation in line with the renormalization group (RG). The main requirement for these functions is to possess the correct UV behavior. This can be achieved by the following ansatz in case of Z 4 : where the momentump is defined asp 2 = (p 2 + q 2 + r 2 + s 2 )/2. The factor 2 appears, because in the integralp approaches for high loop momenta the loop momentum itself. The exponents α 4g and β 4g are determined from the anomalous dimensions of the ghost propagator (δ = −9/44), the gluon propagator (γ = −13/22) and the four-gluon vertex (γ 4g = 2δ − γ = 2/11) from the requirement  Figure 2. Propagator dressing functions for decoupling (straight line) [53] and scaling (dashed line) [78] in comparison to lattice data [79] with β = 6 and lattice sizes of L = 32 (green) and L = 48 (red). α 4g δ + β 4g γ = γ 4g . As second condition serves the IR finiteness of D 4g RG (p, q, r, s) [53]. Solving for the exponents yields For Z 1 the replacement is given by [53]

Input
The four-gluon vertex DSE is calculated with input from other calculations whose results are in good agreement with lattice data. However, this was achieved by the dependence of the corresponding calculations on higher Green functions which were modeled based on the existing information from several sources. The propagators we use for the decoupling case stem from Ref. [53], where the ghost-gluon vertex was dynamically included and an optimized effective three-gluon vertex was used. The results for the propagators, obtained for α(µ) = g 2 /4π = 1, are shown in Fig. 2. Although the ghost-gluon vertex was also calculated there, we employ the bare vertex here. The reason is that it enters only in a static diagram. From the three-gluon vertex, where this is also the case, we know that the influence on the results is minor. We illustrate in Fig. 4 that this also holds true for the four-gluon vertex. On the other hand, for the propagators this is different and the mid-momentum regime of the gluon propagator is affected by this choice [80]. For scaling we use as input data obtained along the lines of [78]. The corresponding dressing functions are shown in Fig. 2.
The ghost-gluon vertex is described completely by one dressing function alone due to the transversality of the Landau gauge: Γ abc µ (k; p, q) := i g f abc P µν (k)p ν D Acc (k; p, q). The gluon momentum is denoted by k and the (anti-)ghost momentum by (q) p and P µν (k) is the transverse projector. Except where mentioned in Sec. 4.1 we use D Acc (k; p, q) = 1.
The three-gluon vertex, on the other hand, has four transverse tensors. However, here we consider only the tree-level tensor and denote the full three-gluon vertex by where the tree-level tensor is given by 3) The dressing function D 3g (p, q, k) contains the non-perturbative information. Note that the restriction to the tree-level tensor is a very good approximation as demonstrated by an explicit calculation of the other dressing functions [58] which were found to be very small. The data we use for the decoupling three-gluon vertex was calculated in Ref. [57] with the propagators discussed before. Again, good agreement with lattice results was found as shown in Fig. 3. For that calculation a model for the four-gluon vertex was required, see Eq. (4.1) below. However, we want to emphasize that coupling this vertex back into the gluon propagator reduces the agreement with lattice results again. Thus we expect that two-loop diagrams are important in the gluon propagator DSE, see also [59][60][61]. The input we have here, on the other hand, can be interpreted essentially as equivalent to existing lattice results. For the scaling solution we calculated the three-gluon vertex from the propagators using the same four-gluon vertex model. Note that although the IR behavior of the corresponding diagram is not correct then, the IR behavior of the three-gluon vertex itself is because it is determined by the ghost triangle diagram. The scaling results for the three-gluon vertex are shown in Fig. 3.
All input quantities are rescaled in the plots, because the shown lattice data is not renormalized. In our calculations, on the other hand, we used the renormalized input data. The importance of consistently renormalized input data is discussed in Sec. 4.  Table 2. Renormalization constants and cutoffs for the used decoupling [53] and scaling input [78].

Ghost box
Since it is expected that the ghost box yields the IR leading contribution to the four-gluon vertex both for scaling and decoupling solutions [33,56,72] we start by a dedicated analysis of this diagram. No iteration is necessary and we can calculate specific kinematic configurations with a very high precision. The ghost box contribution for the three special configurations defined in Tab. 1 is shown in Fig. 4. To obtain the symmetrized results for configuration C, we calculate the (two) distinct permutations and then take the (weighted) average. The other configurations need not be symmetrized since no distinct permutations exist. The symmetrization for configuration C is especially important if the transversely projected decoupling ghost box is calculated. In this case, contributions of different permutations diverge logarithmically in the infrared (IR) but the (weighted) sum approaches a finite value. Thus we confirm the finiteness of the tree-level dressing function beyond configuration A, for which this was already found in Ref. [66]. Indeed our calculations show that the ghost box is IR finite for all momentum configurations.
In Fig. 4 we also illustrate the effect of the transverse projection by showing results obtained from the projector given in Eq. (2.27) without transverse projections. As it turns out, this can have a sizable quantitative but not qualitative effect: The form of the curves stays the same, but the transverse projection can increase (config. B) as well as decrease (config. A and C) the contribution of the ghost box.
The influence of a dressed ghost-gluon vertex was also studied and is shown in Fig. 4. For this we used the ghost-gluon vertex from Ref. [53]. Note that with a dressed ghost-gluon vertex the symmetrization requires the calculation of more diagrams than in the case of a bare vertex. A dressed ghost-gluon vertex always leads to an increase in the region below 1 GeV, but compared to the contributions of the other diagrams, discussed below, this is only a minor effect. It can be shown analytically that the ghost box vanishes completely under transverse projection for configuration A. It was this configuration that was used in Ref. [65] to determine the IR scaling fixed point of the running coupling. We will come back to this in Sec. 4.4 It is noteworthy that the angular dependence, viz., the dependence on the configuration, is stronger if the ghost box is transversely projected.

Full calculation
We will now turn to the results from the self-consistent solution of the truncated four-gluon vertex DSE. For this we take into account the complete momentum dependence and solve the equation by a fixed point iteration as explained in Appendix B.
The results for the four-gluon vertex dressing function are shown in Figs. 5 and 6. The former shows the kinematic configurations A, B and C. The shaded area indicates the angle dependence. To determine it we used the configuration and took the minimal and maximal values for the dressing when varying the angles. However, this area has to be interpreted with a grain of salt, because what we plot as the boundaries of this area is determined by a few extreme points whereas the majority of the points lies around the straight lines. To illustrate this we show three-dimensional plots in Fig. 6 where configurations with the largest angle dependence are shown. Note that the typical angle dependence is much smaller. The main origin of the angle dependence is the gluon box. To check that this is not a numeric artifact, we calculated the gluon Decoupling Scaling Figure 6. Top: angular dependence of the dressing function. The squared momenta and the angle ψr are kept constant: P 2 = Q 2 = R 2 = 160 GeV 2 , ψr = 0. The shown configuration corresponds to the point with the largest angle dependence we found, see also Fig. 5. Bottom: momentum dependence of the dressing function: P 2 = Q 2 = p 2 and R 2 = q 2 . For the angles we chose θq = θr = ψr = π/2 (decoupling); and θq = π/4, θr = π/2 and ψr = 0 (scaling). Note that the lower plots are shown from different viewpoints.
box for momenta and two angles fixed while varying the third angle. This calculation was repeated with increased numeric precision. The results, shown in Fig. 7, clearly illustrate that increasing the numeric precision has no impact. Fig. 7 also shows the effect of a dressed three-gluon vertex compared to a bare one: It strongly enhances the angle dependence. Furthermore, the importance of the symmetrization can be seen. The plot shows the contribution of a single gluon box. If it were already symmetric, the results would be the same independent of varying θ q or θ r . Note that even in the case of bare three-gluon vertices there is a difference, because a single gluon box is not symmetric under the exchange of q and r, see Fig. 14, but only the sum of the three diagrams appearing in the DSE.
The contributions of the individual diagrams to the dressing function are plotted in Fig. 8 for different configurations. As can be seen there, the scaling and the decoupling solution are very similar above 1 GeV. Below 1 GeV, the ghost box starts to dominate the scaling solution due to the IR  divergence, D 4g ∝ p 2 −4κ . For configuration A, for which the ghost box vanishes, the dynamic diagrams become large in the IR but remain finite. For all other configurations, all other diagrams of the scaling solution are insignificant in the IR due to the dominance of the ghost box. For the decoupling solution, since for this tensor there is no IR divergent contribution from the ghost box, all diagrams contribute with a finite value. The mid-momentum regime is dominated by the gluonic diagrams. Interestingly, both triangle diagrams, dynamic and static, are very similar. 3 The two diagrams differ only by the dressings of the three-and four-gluon vertices, see Fig. 1.
So far, the gluonic four-point interactions had to be modelled whenever they were not neglected. In Ref. [60], e.g., a specific model for the four-gluon vertex was used to incorporate the sunset diagram of the gluon propagator DSE. In recent calculations of the three-gluon vertex [57,58], the four-gluon vertex also had to be modelled. In both studies it was found that the four-gluon vertex must be of a certain strength so that the three-gluon vertex DSE converges within the applied truncation scheme. The model employed in Ref. [57] is given by D 4g, dec model (p, q, r, s) = a tanh b/p 2 + 1 D 4g RG (p, q, r, s) (4.1) with D 4g RG (p, q, r, s) defined in Eq. (2.32). It is interesting to see that such a simple form can indeed describe the four-gluon vertex tree-level dressing quite well as we tested by fitting configurations A, B and C. Given that the angle dependence is predominantly weak, the fit for configuration C, shown in Fig. 5, can serve as a good first approximation for the four-gluon vertex in other calculations. The values for the parameters are a = 1.15 and b = 0.63 GeV 2 . However, we emphasize that this function only gives a qualitative representation of the four-gluon vertex. In particular, it describes only the tree-level tensor.

Other tensors
Given the complexity of the four-gluon vertex DSE with its many tensor structures, the calculation of all dressing functions constitutes a further challenge which will not be entered here fully. However, as a first step we can take our solution approximated by the tree-level tensor and calculate other dressing functions. For the three-gluon vertex all transverse dressing functions were calculated in [58] with the result that the tree-level tensor yields indeed by far the most important contribution with the other three tensors at least an order of magnitude smaller. When it comes to the four-gluon vertex, the situation is similar, but with an additional twist: While we find that the other tensor structures we probed are suppressed as compared to the tree-level one over a wide momentum regime, they possess a logarithmic IR divergence.
To assess the size of other dressing functions, we use the results for the tree-level dressing from the full (decoupling) calculation and calculate several other projections. This is no longer a self-consistent solution but should give us at least an idea about the magnitude of such contributions. We consider two classes of other tensors: One that also contains only the metric and no momenta and thus contains the tree-level tensor, and one whose tensors are constructed from momenta only.
The two tensors chosen from the first class are based on Refs. [64,65] but restricted to their transverse parts. They are constructed from a subset of Bose symmetric tensors by orthogonalization. The four-gluon vertex is then written as The basis tensors V abcd i,µνρσ (p, q, r, s) are given in Eq. (C.6). V abcd 1,µνρσ (p, q, r, s) corresponds to the transversely projected tree-level tensor.
Due to the orthogonality of the tensors V abcd i,µνρσ (p, q, r, s), the corresponding dressing functions can be extracted from the four-gluon vertex DSE by replacing the tree-level tensor in the projector (2.27) by the corresponding V abcd i,µνρσ (p, q, r, s). However, in practice we found it easier to project with the non-orthogonalized tensors given in Eq. (C.2), because they are shorter, and calculate the dressings D 4g,Vi (p, q, r, s) from the results.
The two additional dressings D 4g,V2 (p, q, r, s) and D 4g,V3 (p, q, r, s) are shown in Figs. 9 and 10. Again we find no large dependence on the chosen configuration. Most strikingly the magnitude of the two dressings is very small compared to the tree-level dressing, but it becomes larger in the IR where they diverge logarithmically. This divergence, which is due to the ghost box, see Figs. 9 and 10, was also seen in Ref. [66] for a related dressing D 4g,G that belongs to the tensor G abcd µνρσ given by G abcd µνρσ = (δ ab δ cd + δ ac δ bd + δ ad δ bc ) × (δ µν δ ρσ + δ µρ δ νσ + δ µσ δ νρ ). G can be written as a linear combination of the tensors V 2 and V 3 , see Eq. (C.2). Since the results for D 4g,V2 and D 4g,V3 are very similar, D 4g,G resembles the two as well. In particular we do not find an enhancement in the mid-momentum regime as it was found in Ref. [66]. As we checked explicitly, the source of this enhancement lies in the truncation scheme employed in Ref. [66], where all renormalization constants in front of the loop diagrams were set to 1. One can see in Figs. 9 and 10 that the contributions of the individual diagrams are by no means small and the resulting dressing is only small in the mid-momentum regime because of delicate cancellations. Since all diagrams contain either Z 1 or Z 4 , which have very different values, see Tab. 2, the renormalization constants naturally play an important role in this and should be taken into account properly. As a second example we consider a tensor constructed from momenta only. This case serves to investigate if tensors from another class as before have a different behavior. As a Bose symmetric representative of such a class we choose the tensor P abcd µνρσ (p, q, r, s) given in Eq. (C.6). It is constructed from P abcd µνρσ = (δ ab δ cd + δ ac δ bd + δ ad δ bc ) × s µ r ν q ρ p σ + r µ s ν p ρ q σ + q µ p ν s ρ r σ p 2 q 2 r 2 p 2 (4.4) , config. C ghost box gluon box gluon triangle swordfish dynamic triangle Figure 9. Results for D 4g,V 2 (left) and individual contributions to configuration C (right).   by transverse projection and normalization to the tree-level, see Appendix C. As it turns out the behavior of the corresponding dressing is indeed different from that of D 4g,V2 and D 4g,V3 . In particular the angle dependence is more pronounced as can be seen in Fig. 11. Besides the height of the bump in the mid-momentum regime also the sign depends on the configuration. We present no results for configuration A, since the tensor P is then purely longitudinal.
Finally we compare the contributions of different dressing functions in Fig. 12. Clearly, the tree-level dressing function is dominant. However, this is due to the tree-level diagram itself, which contributes with the constant value Z 4 , see Fig. 8. The loop diagrams only account for comparatively small changes around the value of the renormalization constant Z 4 . Taking this into account, the contribution of the loop diagrams are similar for all considered dressings. One exception is that only the non-tree-level dressing functions diverge logarithmically in the IR. However, this divergence sets in at very low momenta at around 100 MeV and the dominant dressing over a wide range of momenta is that of the tree-level tensor.

Running coupling
From the four-gluon vertex a renormalization group invariant running coupling can be defined, as it can be from any other vertex [22,72]. Up to now the couplings derived from the ghost-gluon, the three-gluon and the four-gluon vertices were calculated, see, e.g., [32,53,57,58,65] for results from DSEs. They are given by [72] α MM (p 2 ) = α(µ 2 ) G 2 (p 2 )Z(p 2 ), (4.5a) where α(µ 2 ) = g 2 /4π is the running coupling at the renormalization scale µ. For the arguments of the three-and four-gluon vertices a generic scale p 2 was given. Which kinematic configuration is chosen is in principle free. For the three-gluon vertex we choose the symmetric point and for the four-gluon vertex configuration C. Note that the denominators of the three-and four-gluon vertex couplings are not unity, because we work here in the MiniMOM scheme where G 2 (µ 2 )Z(µ 2 ) = 1. Choosing, e.g., Γ 3g (µ 2 ) 2 Z 3 (µ 2 ) = 1 would correspond to a different scheme. Hence we have to take this factor into account.
In Fig. 13 we show the different running couplings. For large momenta they all agree as they should with only small deviations. At low momenta the four-gluon vertex running coupling vanishes like p 4 . In the scaling case all couplings exhibit an IR fixed point [72]. For the four-gluon vertex this was confirmed in Ref. [65]. However, configuration A was used which, as discussed in Sec. 4.1, vanishes upon transverse projection. Thus this configuration is not suited for calculating the running coupling. Since the existence of the IR fixed point does not depend on the angular configuration, we believe that is a shortcoming of the truncation. Indeed there are two IR leading diagrams we neglected that may yield a contribution that does not vanish in the IR and provides a non-vanishing value for the running coupling. For configuration C we extract a value of α 4g (0) = 0.000 42. This is in accordance with the findings of Ref. [65] that the fixed point value of the four-gluon vertex is much lower than that of the ghost-gluon vertex, which is α MM (0) ≈ 2.97 [25]. For the symmetric point of the three-gluon vertex we extract α 3g (0) = 0.0032. Our value for α 3g (0) deviates from the value found in Ref. [58], α 3g (0) 0.0016, but this is due to different ghost propagator input.

Summary and conclusions
In the past the two-and three-point functions of Yang-Mills theory were intensively scrutinized with functional equations. However, within the employed truncations these calculations still relied on model input for higher Green functions. Using lattice data, these models could be tuned such that neglected contributions could be effectively taken into account which led to good agreement with lattice results. Thus, as far as two-and three-point functions are concerned, we have quantitatively reliable input for other calculations such as those performed here.
Employing the by now common truncation to the UV leading diagrams, which also contain IR leading diagrams, the DSE of the four-gluon vertex does not rely on any model. This is, within this truncation scheme, a unique feature among the primitively divergent Green functions. However, given the tensorial complexity of the vertex, we solved the DSE self-consistently only for the tree-level tensor. A posteriori we confirmed for some additional dressing functions that their magnitude is much smaller than that of the tree-level tensor over a wide momentum regime. In the IR, on the other hand, the tree-level dressing is in the decoupling case constant, whereas other dressings can diverge logarithmically. These divergences set in at about 100 MeV. Since a four-gluon vertex is within a functional equation for 1PI functions contracted with at least two gluon dressing functions, which are IR suppressed, these divergences most likely do not have a strong effect on other Green functions as long as no other IR divergent expressions appear in the integrand.
For the dominant tree-level structure we found that the box diagrams are almost negligible and the triangle and swordfish diagrams yield the largest contributions. However, for other tensors this is no longer the case. Especially the gluon box can have a sizable impact and the ghost box leads to the IR divergence mentioned before.
It is important to note that the fact that other dressings functions are so small in the midmomentum regime is in no way trivial, as it comes from cancellations between all the diagrams taken into account. If it turns out that also dressing functions beyond those we investigated here follow this pattern, the four-gluon vertex can be well approximated by one tensor only. This would alleviate its use in future studies, like its effect in the gluon propagator or three-gluon vertex DSEs, considerably. As a first approximation we provide a fit that describes the dressing qualitatively well.
An important aspect of our calculations was that we used only the transverse subspace, since it is sufficient for the Landau gauge. By studying the ghost box diagram explicitly also without transverse projection we found that there can be a considerable influence of this restriction which most likely also exists for other diagrams. Thus in future studies this restriction should always be taken into account.
From the four-gluon vertex a running coupling can be extracted. Qualitatively it behaves as expected, viz., it agrees in the perturbative regime with the couplings from the ghost-gluon and threegluon vertices and then turns towards zero in the IR for decoupling and towards an IR fixed point for scaling. We confirmed that the value of this fixed point is very small compared to that of the MiniMOM coupling.
The four-gluon vertex was the last primitively divergent Green function of Landau gauge Yang-Mills theory for which a self-consistent solution was lacking. Its calculation constitutes an important step towards a fully self-contained description of Yang-Mills Green functions from functional equations and will enable the study of its impact on the propagators and the three-gluon vertex.

A Color calculations
To calculate color traces, we employ the following two well-known identities: If the trace contains six structure constants, it is sometimes necessary to insert the Jacobi identity We show now that C 7 , defined in Eq. (2.22), is not linearly independent of the tensor set given in Eq. (2.21) for N c < 4. To start, consider the set of tensors C 1 , . . . , C 5 and C 7 . The scalar product C i , C j of two tensors is given by the trace over the color indices. Thus, the metric tensor of a six-dimensional color space defined by C 1 , . . . , C 5 and C 7 is given by where C i is defined as the normalized C i , viz., .

(A.4)
Hence, if N c = 3 or N c = 2, det(g color 6 ) = 0. If the determinant of the metric tensor is zero, the tensors (i.e., C 1 . . . C 5 and C 7 ) are linearly dependent. One can similarly show that the determinant of the metric tensor of the color space spanned by the tensors C 1 . . . C 5 is not zero, det(g color 5 ) = 0 . Thus, C 7 can be expressed in terms of C 1 . . . C 5 : (

B Technical details of the calculation
The four-gluon vertex was calculated with the CrasyDSE framework [69]. In a first step, we calculate all static diagrams, the diagrams that do not depend on the four-gluon vertex (the diagrams in the second line in Fig. 1). Calculating on cores of Intel Xeon E5-2670 processors, the first step takes typically 30 000 core hours. For the actual iteration process we only need to calculate the swordfish and the dynamic triangle diagrams, where one iteration step typically takes 10 000 core hours. Fortunately, the four-gluon vertex DSE converges relatively fast within 6 to 8 iteration steps. We perform the integration as detailed below. The integration itself is done by a standard Gauß-Legendre quadrature. Unfortunately we cannot integrate out any variables analytically and have to perform all four integrations numerically. We use spherical coordinates given by  The integral measure is then (with K = |k|) Exploiting that all scalar products of k and s, r or q, and thus the kernels, depend on cos(φ k ) only, the φ k integration can be simplified: 2π 0 dφ k · · · = 2 π 0 dφ k . . . . We split each integration into several integration regions in order to not integrate over the singularities arising from the internal propagators. We choose the momentum routing such that the loop momenta of the gluon box and the ghost box are given by The routing is illustrated in Fig. 14. The integrand diverges if k 2 1 = 0, k 2 2 = 0, k 2 3 = 0 or k 2 4 = 0. Due to the choice of coordinates, this leads to four integration regions in K 2 , three in θ k , two in ψ k and one in φ k . The momentum routing of the swordfish and the triangle diagrams can be chosen such that only loop momenta given in Eq. (B.3) appear. This allows us to use the same quadrature for all diagrams. For the static diagrams we use 30 (12) integration nodes for every momentum (angular) integration region (for the dynamic diagrams we use up to 15 integration nodes per angular integration region). Thus, for each of the 15 3 · 7 3 = 1 157 625 external grid points, we need to evaluate 1 244 160 internal grid points. The most complicated kernel is that of the gluon box. To evaluate the gluon box kernel once, we need about 50 000 multiplications even though we optimized the kernel with a specialized Mathematica [81] algorithm.
In the evaluation of the kernels values for the variables of the four-gluon vertex may appear that are outside of the grid. In this case an appropriate extrapolation must be performed. For the angles (θ q , θ r , ψ r ), it can happen that they are (slightly) higher than the highest angular grid point, since we did set the highest grid point slightly below π to avoid the appearance of a (well-defined) 0/0. In that case we simply approximate the dressing function by the boundary value.
The squared momenta, on the other hand, can be below as well as above the range defined. In the former case, the IR extrapolation, we again use the boundary value. This is reasonable if the IR behavior of the dressing function is insignificant for the iteration, or the dressing function is constant in the IR. The latter is approximately the case in the decoupling solution (see, e.g., Fig. 8), the former in the scaling solution: The ghost box is the leading diagram of the scaling solution in the IR. Since the ghost box does not depend on the four-gluon vertex, it does also not depend on the extrapolation. Thus the IR region is not affected by the extrapolation of the four-gluon vertex.

(B.4)
We find that the highest calculated UV points show the expected UV behavior given by Eq. (B.4).
In addition, Eq. (B.4) connects the extrapolated momentum region with the non-extrapolated region smoothly. Therefore, the UV extrapolation by Eq. (B.4) is a reasonable extrapolation. However, in some isolated cases it can happen that the results deviate. One such example is configuration A with p 2 equal to the highest point on the momentum grid. For specific permutations the angle configuration corresponds to one of the configurations with a large angle dependence, see Fig. 6. Thus the extrapolation gives a value that is off. However, such configurations are rare and we expect them to have a negligible influence on the total calculation. An extrapolation is also necessary for some points in the symmetrization process. However, the UV extrapolation function is valid for the dressing function D 4g (p, q, r, s) and not L(p, q, r, s) given in Eq. (2.27). Thus we work with D 4g (p, q, r, s) = Z 4 + 3L(p, q, r, s) (B.5) which should approximate the dressing function well enough to allow a good extrapolation. The final dressing function is obtained from averaging over all 12 D 4g .
Here we are only interested in the transverse part of the four-gluon vertex and thus construct a basis from the transversely projected tensors. Combining this with a Gram-Schmidt orthogonalization we obtain the following basis, where all indices and arguments are suppressed and T denotes the transverse projection of all four Lorentz indices: As tensor of the second class we consider P abcd µνρσ = (δ ab δ cd + δ ac δ bd + δ ad δ bc ) × s µ r ν q ρ p σ + r µ s ν p ρ q σ + q µ p ν s ρ r σ p 2 q 2 r 2 p 2 , (C.4) which is also Bose symmetric. The tensor P is orthogonal to V 1 due to the color part but not to V 2 and V 3 . For estimating the magnitude of the corresponding dressing functions orthogonality to the tree-level tensor suffices. Hence, we simply use the transversely projected tensor (suppressing indices again) P = T P . (C.5) For meaningful comparisons the tensors need to be normalized. Thus, we introduce normalization factors: V i = N i V i , P = N P P (C.6) To simplify notation, we denote the norm of a tensor by X = X abcd µνρσ X abcd µνρσ .
A straightforward choice would be to normalize the tensors to 1, viz., V i = 1 and P = 1. However, then V 1 does not coincide with the transversely projected tree-level tensor. Thus we use N 1 = 1. In order to make the dressing functions comparable, we demand which leads to , and N P = V 1 P . (C.8) The signs of the tensors are in principle arbitrary and we chose them such that the corresponding dressing functions have the same sign as the tree-level. Note that via the normalization factor each tensor gets a factor of g 2 . In general, the normalization factors depend on the momenta. However, for a fixed momentum configuration with one momentum scale, as used in our plots, the norms are constant for purely dimensional reasons. Thus, the normalization factors are also constant then. Finally we mention the influence of this choice of normalization on Fig. 12. To obtain tensors normalized to 1, we need to divide the tensors V 1 , V 2 , V 3 and P by V 1 . Since the V 1 is constant for a specific configuration, the relative values of the dressing functions in Fig. 12 were identical if we had chosen to work with tensors normalized to 1. To be precise, a plot for tensors normalized to one can be obtained by multiplying all values of the dressing functions by V 1 = 7 2 g 2 N c 3(N 2 c − 1) ≈ 51.44 g 2 for configuration C.