Correlation functions of Landau gauge Yang-Mills theory

Correlation functions of Yang-Mills theory in the Landau gauge are calculated from their equations of motion. The employed setup is completely parameter-free and leads, within errors, to good quantitative agreement with corresponding lattice results for the ghost and gluon propagators as well as the ghost-gluon and three-gluon vertices. Also the four-gluon vertex is calculated. The present setup allows for the first time for a unique subtraction of quadratic divergences in the gluon propagator Dyson-Schwinger equation. Thus, there are no ambiguities which can arise due to the use of models or auxiliary workarounds. In addition, several self-tests of the results are described that allow assessing the truncation error in a self-consistent way. This enables a new perspective on how to identify limitations of the present setup and develop future improvements.


I. INTRODUCTION
Quantum chromodynamics (QCD) describes the strongly interacting part of the Standard Model and exhibits a rich phenomenology. Describing it from first principles is a challenge for contemporary physics. Functional methods are one possible approach for this. They are formulated in terms of the correlation functions of quarks and gluons and constitute a system of nonlinear equations. To solve them, some form of approximation is required, and often models are used to fill the gaps. This has led to a successful description of many bound states and allowed exploring the phase structure of QCD, see [1,2] and [3,4] for reviews, respectively. In recent years, though, some of these gaps have been closed, and a coherent picture is emerging of how functional equations can lead to a self-contained description without any modeling. Thus, it seems within reach to establish a direct link between the QCD action and its phenomenology with functional equations. The aim of this work is to combine several of the recent individual advances using equations of motion and provide a state-of-the-art solution for the correlation functions of Yang-Mills theory. Understanding Yang-Mills theory, which describes only the gluonic part of QCD, is a cornerstone for the treatment of full QCD, but solving the corresponding functional equations faces a few challenges which were overcome only recently. Some of these challenges were of a technical nature which required to adopt and develop the correspond tools. Others are related to problems which can be alleviated by adjustments of the employed models or other modifications of the equations. However, for a self-contained calculation this freedom is lost, as by definition no free parameters are allowed. Consequently, these problems have to be * markus.huber@physik.jlug.de handled properly. Functional equations are used with various gauges, but here solely the Landau gauge is considered. It has many advantageous features which make it a very convenient choice. Propagators and vertices were calculated with different types of functional equations using a variety of different approximations. Other methods, like lattice simulations or effective models, have been used as well. However, to test the robustness and reliability of the presented results, not only comparisons with other methods will be performed, but also some self-tests. Such tests are crucial to make the method independent from the availability of other results. In addition to elucidating structures and mechanisms qualitatively, the method then becomes also quantitatively predictive starting from first principles. The final test of results for gauge-dependent quantities consists in calculating a gauge-invariant quantity. When, as is the case here, no free parameters exist, the selfconsistency of the whole setup is decisive, as any inconsistencies cannot be covered up anymore by tuning a model. For Yang-Mills theory a natural, though not easy example for such a gauge-invariant quantity is glueballs. While fore real glueballs the effects of quarks still need to be added, we can for pure Yang-Mills theory settle with a comparison to corresponding lattice calculations. The results presented here were used to this end in Ref. [5] where the masses of scalar and pseudoscalar glueballs were calculated. The good quantitative agreement obtained is an important indication that the results from the employed truncation have reached the necessary quantitative precision, at least for this task. The encouraging agreement with other methods and the successful calculation of gauge-invariant quantities should not hide the fact that there is still room for improvement. Having reached the present level of accuracy, one should try to identify the remaining shortcomings as a next step. They might be small, but depending on the specific task, they might still play a relevant role. Thus, it is important for the future improvement of truncations to understand their sources and how to get rid of them. In this respect it should be stressed that qualitatively this process is now different than in smaller truncations because there are no ambiguities left. A prime example of this are quadratic divergences which are discussed in detail in Sec. IV. In the past, their removal in the gluon propagator equation of motion had a quantitative influence on the result which made it difficult to disentangle the effects of different parts of the truncation. Hence, using the present truncation as a starting point, the impact of various improvements can be studied in a much clearer way. To this end, detailed comparisons and analyses are presented in Sec. V and Sec. VI. The remainder of the article is organized as follows. Yang-Mills theory in the Landau gauge and the correlation functions calculated in this work are introduced in Sec. II. The employed functional equations are discussed in Sec. III including details on the truncation. Their renormalization is explained in Sec. IV. Results are presented in Sec. V and compared to results from other methods in Sec. VI. Sec. VII contains a discussion of the results and Sec. VIII the conclusions. Results for the system using an alternative equation for the ghost-gluon vertex, the renormalization procedure for the scaling case and numerical details are deferred to appendices.

II. CORRELATION FUNCTIONS OF YANG-MILLS THEORY
The Lagrangian density of Yang-Mills theory with the gauge group SU (N ) and the gauge fixed to linear covariant gauges is 1 where the field strength tensor and the covariant derivative are given by D rs µ = δ rs ∂ µ + g f rst A t µ .
respectively. The gluon field A and the ghost (anti-ghost) field c (c) live in the adjoint representation. The fundamental generators obey the relations Tr{T r T s } = 1 2 δ rs .

A. Propagators
The gluon propagator is given by D ab µν (p) = δ ab D µν (p) = δ ab (D T µν (p) + D L µν (p)), 1 The conventions in this article follow those of Ref. [6]. D T µν (p) = g µν − p µ p ν p 2 D(p 2 ), D(p 2 ) = Z(p 2 ) p 2 , (7) In the Landau gauge, the propagator is completely transverse, viz., ξ = 0. The gluon propagator is then completely described by the scalar part D(p 2 ) = Z(p 2 )/p 2 . The propagator of the ghost field is B. Three-point functions Landau gauge Yang-Mills theory has two three-point functions, the ghost-gluon and the three-gluon vertex. The former has two dressing functions, of which in Landau gauge only one is relevant. Thus, the vertex is written as Γ Acc,abc,T µ (p 2 ; p 1 , p 2 ) = i g f abc D Acc (p 2 2 ; p 2 1 , p 2 2 )P µν (p 2 )p 1ν .
Note that although the transverse projector is put explicitly, this would not be necessary since the vertex is always contracted with a transverse projector anyway. The fields in the superscript indicate the order of the arguments. Here and for other vertices, all momenta are incoming. The bare vertex is Γ (0),Acc,abc µ (p 2 ; p 1 , p 2 ) = i g f abc p 1ν .
The bare vertex is given by In the following, only a dressing function C AAA (p 2 1 , p 2 2 , p 2 3 ) of the tree-level tensor is considered, viz., for the full vertex is used. The color structure employed for these vertices is fixed to the antisymmetric structure constant f abc . In principle, SU (N ) also has the symmetric structure constant d abc for ; see [7] and [8], respectively, for more details. Kinematic points of interest for three-point functions are the symmetric point p 2 1 = p 2 2 = p 2 3 , the lines of two equal momentum squares p 2 i = p 2 j , the lines of orthogonal momenta pi · pj = 0, and the points of one vanishing momentum p 2 i = 0. Similar points exist for four-point functions with k = p1 + p2, p = p2 + p3 and q = p3 + p1.
N > 2. However, due to the charge invariance of QCD, only the totally antisymmetric f abc is allowed [9,10]. Three-point functions depend on two independent momenta from which three kinematic variables can be constructed. A typical choice are two momentum squares and an angle. However, to exploit potential symmetry properties of the vertices and for reasons of technical simplifications, another set motivated by the symmetry properties of the S 3 permutation group is used here [7]. This basis makes the Bose symmetry of the three-gluon vertex manifest, but it is also useful for the ghost-gluon vertex which lacks this symmetry. From the three momenta p 1 , p 2 , p 3 , one constructs the following variables [7] which are a singlet (S 0 ) and a doublet (a, s) under the permutation group S 3 : S 0 can take all positive values and a and s are restricted to a unit disk, a 2 + s 2 ≤ 1, see the left plot of Fig. 1. Computationally, it is thus advantageous to use radial coordinates for a and s: For clarity, the arguments of the dressing functions will be given in the following by the momenta and not the actually used kinematic variables. The inverse transformation is C. Four-point functions Finally, the four-gluon vertex has to be discussed. It has 136 tensors in Lorentz space, 41 of which are transverse [8]. In color space, there are nine independent tensors for general SU (N ), but this number reduces to eight for SU (3) due to an additional identity [11] and to three for SU (2) due to the absence of the symmetric structure constant. However, for SU (3) only five color tensors are required since the other three decouple [12,13]. The tree-level four-gluon vertex reads For the full vertex only a dressing function F AAAA (p 1 , p 2 , p 3 , p 4 ) of the tree-level tensor is taken into account: The dressing function F AAAA (p 1 , p 2 , p 3 , p 4 ) depends on three independent momenta out of which six variables can be formed. They can be organized according to the S 4 permutation group into a singlet, a doublet and a triplet [8]. From the three independent momenta p 1 , p 2 and p 3 , the following momentum combinations are constructed: The momentum squares k 2 , p 2 and q 2 are used together with the scalar products between different momenta, to define the singlet variable 2 the doublet variables and the triplet variables The doublet variables are similar to the case of threepoint functions, except that a and s are restricted to the interior of a triangle [8], see the right plot of Fig. 1. Also in this case, a and s are expressed via radial coordinates: For symmetry reasons, the angle η is measured from the negative s-axis and not the a-axis as for the three-point functions. The different shape of the domain of a and s is encoded in an angular dependence of the radial coordinate ρ(η). It is convenient to rescale it to the interval [0, 1] asρ where ρ max is given by · denotes the floor function. The dependence of the four-gluon vertex on six kinematic variables is slowing down calculations considerably. As an approximation, only the singlet and doublet are taken into account here by setting ω 1 = ω 2 = ω 3 = 0. This entails p 2 1 = p 2 2 = p 2 3 = p 2 4 = S 0 . A nontrivial angular dependence enters via the doublet variables: Even with this approximation, the calculation of the fourgluon vertex dominates the computing time.

III. EQUATIONS OF MOTION
In this section, the equations of motion from the 1PI effective action (Dyson-Schwinger equations) and higher nPI effective actions, which will be used to calculate propagators and vertices, are discussed. For simplicity, we refer with nPI effective action always to higher nPI effective actions with n > 1 and denote their equations of motion as EOMs, while the ones from the 1PI effective action are denoted by DSEs. The equations of both cases are similar and can be treated numerically with the same methods. A difference consists in the structures of the complete systems of equations: DSEs form an infinite set of equations of finite size, while nPI effective actions lead to a finite number of equations with (possibly) infinitely many terms. Such differences become important when considering how to truncate the equations. For more details beyond the short overview presented here, see, e.g., [6,[14][15][16][17][18]. After a general discussion of the equations, the specific equations for the two-, three-and four-point functions are presented and how they are truncated.

A. Dyson-Schwinger equations
The 1PI effective action Γ[Φ] depends on the classical fields only. In the shorthand notation employed here, Φ is the classical super-field that represents all fields, i.e., Φ = {A,c, c}. It is related to the underlying action S[φ] via a Legendre transformation of the generating functional of connected correlation functions W [J], depending on the corresponding sources J of the quantum fields φ, with Φ = δW/δJ = φ J : Summation as well as integration over repeated indices, which represent field-types and all internal and space/momentum variables, are understood. DSEs can be derived from the invariance of the path integral under translations of the fields which lead to (see, e.g., [6,14,15,17] for details) The field-dependent (as indicated by the index J) propagator D ij J is given by The physical propagators are obtained by setting the sources to zero: D = D J=0 . From the master equation 38, DSEs for all correlation functions can be derived by differentiating with respect to the appropriate fields and setting the sources J to zero. This leads to infinitely many coupled equations. The equation for an n-point function does not only depend on lower correlation functions but also on (n + 1)-and, due to the presence of a four-point function in the action, possibly on (n+2)-point functions. However, each equation has a finite number of terms, as can be directly inferred from Eq. (38) [27] which is also used for solving the equations numerically.

B. Equations of motion from nPI effective actions
The 1PI effective action is formulated as an expansion in the classical field Φ with the fully dressed 1PI correla-tion functions of the theory, denoted by Γ (i) , as expansion coefficients. nPI effective actions differ insofar, as propagators and vertices up to order n are treated on the same footing as the fields and are taken into account via additional Legendre transformations. Vertices with i > n only appear bare. In analogy to the field variable, additional sources R (i) , 2 ≤ i ≤ n, are introduced: The nPI effective action is then obtained from corresponding Legendre transformations and depends in addition to the classic field Φ also on the propagator D and potentially vertices Γ (i) , where arguments have been suppressed on the right-hand side. Differentiating the effective action and setting all sources to zero leads to stationarity conditions which correspond to the EOMs: There are as many equations as there are Legendre transformations with respect to sources J and R (i) . In contradistinction to DSEs, these equations have infinitely many terms. In practice, a loop expansion of the nPI effective action is performed. Truncating it renders the numbers of terms finite. In addition, using the EOMs of higher correlation functions, cancellations in an equation can be realized [16,18] that lead to simplifications. An nPI effective action for Yang-Mills 3 theory can be written as S ij is the field-dependent inverse bare propagator, defined as For the full propagators D ij , their field content was made explicit with the superscripts. Γ is split into parts with (Γ int ) and without (Γ 0 ) bare vertices. Here, we will consider the 3PI effective action at three-loop level which is required for a self-consistent truncation including threepoint functions. The corresponding Γ is given in Fig.  2 The truncation of a system of DSEs specifies how to reduce the infinite system of equations to a finite set of equations. This can be done by setting certain correlation functions to zero or by providing models for them. Also EOMs can be truncated in this way, but nPI effective actions offer a more systematic way via a loop expansion. One then obtains a concrete form for the nPI effective action from which all EOMs are derived. Changing the truncation leads then to consistent changes in all EOMs and manifests the relation between different terms. Here, we will discuss truncations of DSEs and EOMs and describe the system of equations that was used for the calculations. It should be noted that for the self-consistent treatment of an m-point function, an l-loop expansion of the nPI effective action with l ≥ n (and trivially m ≥ n) is needed.
The goal is to calculate all primitively divergent correlation functions of Yang-Mills theory: The gluon and ghost propagators, the ghost-gluon vertex, the three-gluon vertex and the four-gluon vertex. The lowest nonprimitive ones, the two-gluon-two-ghost and four-ghost vertices, were analyzed in Ref. [13]. Within the employed kinematic approximations, their impact on the DSEs of primitively divergent correlation functions was found to be negligible. Thus, they are neglected here as are fiveand higher n-point functions. In the following, the equations for individual correlation functions and details of the employed truncations are discussed. The simplest equation is that for the ghost propagator. Its full DSE is depicted in Fig. 3. For nPI effective actions with n ≥ 3, some diagrams in the corresponding EOM can be summed up such that the EOM equals the DSE [16,18]. In all calculations, the full ghost DSE is used. Also for the gluon propagator resummations can be performed and its EOM looks similar to the DSE. The latter is depicted in Fig. 3. The EOM from an nPI effective action with n ≥ 4 has the same form [18]. For the threeloop truncated 3PI effective action, the only difference is the replacement of the dressed four-gluon vertex in the sunset diagram by a bare one. Although within a certain approximation [28] or for three-dimensional Yang-Mills theory [29] it was found that the sunset is subleading in the gluon propagator DSE, this is not the case here. On the contrary, it is found that a dressed vertex is required for the stability of the solution. It should be stressed, though, that this statement is made within the given truncation. If the truncation were changed, for example, by including all tensors of the three-gluon vertex, this would need to be tested again. For now, the dressed four-gluon vertex is always included in the calculations.
In summary, at the two-point level all equations are untruncated DSEs and thus exact. The complete DSE of the three-gluon vertex, depicted in Fig. 4, contains 14 diagrams. Only five of them are perturbative one-loop diagrams. The sixth one-loop diagram, the ghost swordfish diagram, is perturbatively a two-loop diagram. As mentioned above, it is neglected here due to the negligible impact it has [13]. Before the question of the impact of the two-loop diagrams is discussed, we compare the one-loop truncated DSE to the EOM of the three-loop truncated 3PI effective action which is also depicted in Fig. 4. This equation is very similar to the one-loop truncated DSE (without ghost swordfish) with the difference that all three-/four-point functions are dressed/bare. Obviously, the one-loop truncated DSE is not Bose symmetric, while the EOM is. This can be remedied for the DSE by symmetrizing the result [7,30], see also App. C. However, it is convenient to symmetrize also the result from the EOM to cancel subleading numeric effects arising from having to make a specific choice for the legs with momenta p 1 and p 2 .
Despite the similarities of the equations we will see below that they yield very different results and the EOM result is closer to the result from the DSE with two-loop diagrams included. The calculation of two-loop diagrams of three-point functions is complicated. To assess their importance, they were calculated within a one-momentum approximation.
The diagrams with the two-ghost-two-gluon and fivegluon vertices were not included. From previous calculations using one-loop truncations, it is known that the leading kinematic dependence comes from the momentum scale [7,30,31]. When using the S 3 variables, this is reflected in the dependence on the singlet S 0 and only a small dependence on the doublet variables a and s. In the left plot of Fig. 5, this dependence is shown via the small band. Thus, neglecting the angular dependence of the vertex should provide a reasonable first guess to test the impact of two-loop diagrams and simplifies the calculation. Still, the computational cost increases considerably compared to a one-loop truncation, because seven instead of three integrations have to be performed numerically.
In the right plot of Fig. 5, results for the three-gluon vertex from different equations and truncations are compared: from the two-loop truncation of its DSE described above and the EOM of the three-loop truncated 3PI effective action. Since the former was evaluated using a kinematic approximation, the EOM was evaluated with and without the same approximation to confirm how small the effect of it is. The equations were solved for a fixed input, for which a solution of the full system of equations was used. It would be interesting to compare these solutions to those of the one-loop truncated DSE. However, with the given input, no solution is obtained in that case because the gluon triangle is too strong compared to the swordfish diagrams and the three-gluon vertex equation diverges. This imbalance depends on the details of the input and the situation can be different for variations of it. In this respect, it is also important that no RG improvement terms are included for the bare vertices. Such terms were found to have a substantial impact on the solution by shifting the zero crossing to lower momenta [7,30]. Conversely, without RG improvement terms, the zero crossing is at scales above 1 GeV which makes such solutions very different to the ones shown in Fig. 5. For the EOM setup, solutions with and without the same kinematic approximation that was used for the two-loop DSE are shown, but the difference is marginal. The solution of the two-loop DSE is very similar to that of the EOM, with small differences around 2 GeV and below. Within errors, both results agree with lattice results. Thus, since it is much easier to evaluate and one does not need to resort to kinematic approximations, for the three-gluon vertex the EOM is used instead of the

DSE.
The ghost-gluon vertex has two different DSEs that are distinguished by the leg attached to the bare vertex. For a detailed discussion of the differences see Ref. [6,29]. In Fig. 6, the DSE with the ghost leg attached to the bare vertex and the EOM of the three-loop expanded 3PI effective action are shown. The DSE where the gluon leg is connected to the bare vertex in the diagrams contains two-loop diagrams. It is not only more difficult to calculate, one would also need the two-ghost-two-gluon and the two-ghost-three-gluon vertices which we do not include here. The other DSE consists only of one-loop diagrams but also contains the two-ghost-two-gluon vertex. However, the influence of this diagram was found to be very small at least under certain approximations for the four-point vertex [13]. One can compare this truncation to the EOM of the three-loop truncated 3PI effective action, also shown in Fig. 6. It looks very similar to the truncated DSE with the exception that all vertices are dressed. By comparing results from both equations, one can get an estimate of the truncation error. In the fol- lowing, for consistency with the three-gluon vertex, the EOM is used. The effect of using the DSE instead is discussed in Appendix A. In summary, the vertex itself still depends on the choice of truncation, but the other correlation functions are hardly affected. The DSE for the four-gluon vertex truncated at one-loop level is depicted in Fig. 7. For the calculations, these diagrams without the ghost triangle, the impact of which was found to be small [13], are used. The four external legs lead to a proliferation of combinatoric possibilities of diagrams. However, when solving the equation, each diagram type is calculated only once. From this result the permutated diagrams can be obtained, see Ref. [32] and App. C for details. The EOM from the four-loop expanded 4PI effective action looks very similar to the one-loop part of the DSE but with all vertices dressed and an additional diagram type, a swordfish with two two-ghost-two-gluon vertices.

IV. RENORMALIZATION
Due to their nonperturbative nature, the renormalization of DSEs can be an intricate issue. In particu-lar, when several correlation functions are involved, this problem becomes severe. Also, perturbative resummation is closely related to a proper renormalization. This becomes apparent when tracing the origins of different contributions to the perturbative series of the one-loop resummed expression. To illustrate this, let us consider the perturbative, one-loop resummed expression β 0 is the first coefficient of the beta function and γ the anomalous dimension, e.g., for the gluon. The contribution proportional to g 2 ln p 2 /µ 2 comes directly from the one-loop diagrams. The higher order terms, however, are combinations from different sources. In Ref. [6], it was shown explicitly for φ 3 theory how the terms of order g 4 ln 2 p 2 /µ 2 are a combination of two-loop contributions and one-loop contributions including counter terms. The same happens for QCD and the employed truncation needs to fulfill certain criteria to be able to recover the one-loop resummation. In particular, renormalization needs to be done properly, as the counter-terms contribute directly via the renormalization constants. Thus, a naive two-loop calculation alone is not sufficient, because also the renormalization constants in front of the one-loop diagrams are crucial. A consistent inclusion of two-loop diagrams in the gluon propagator DSE respecting this was realized for the first time in Ref. [13]. The calculation of the renormalization constants will be explained later.
In most previous calculations, the lack of perturbative self-consistency lead to problems with regard to the resummed perturbative behavior. This problem is not specific to Yang-Mills theory but applies to full QCD as well.
Despite variations in the reasoning, all workarounds led to the introduction of artificial terms in the equations. For example, the momentum-independent renormalization constant can be replaced by a momentum-dependent function, e.g., [7,30,[32][33][34][35][36], or vertex models can be modified appropriately to cancel the renormalization constants and correct the anomalous dimensions, e.g., [37][38][39]. In both cases, one gets rid of the cutoff-dependent renormalization constants in a way that the RG behavior of the equation is respected. However, there is always a model-dependent aspect to this which one would like to overcome if one aims at quantitative predictive power. For the Yang-Mills propagators, it was demonstrated in Ref. [13] how this is achieved. Here, it is extended to all other primitively divergent correlation functions. The crucial ingredient is the interplay of renormalization and higher-loop contributions, see Ref. [6] for details.
We start the discussion of renormalization with the simplest equation, the ghost propagator DSE. Its renormalization is realized via subtraction at vanishing momen-tum, where Σ G is the ghost self-energy, Z 3 the ghost wave function renormalization constant and x = p 2 . G(0) can be varied to obtain different solutions [40]. In contrast, in functional renormalization group (FRG) calculations the initial conditions for the gluon propagator can be varied to achieve the same effect [41]. For G(0) → ∞, the so-called scaling solution is obtained [33,[42][43][44][45] for which the dressing functions obey power laws. The corresponding exponents are related to each other [43,44,46,47], unique [48][49][50][51][52] and can be calculated analytically in terms of κ = 0.595 [43,44]. The qualitative behavior, including the scaling relation for the propagators, is in agreement with the Gribov-Zwanziger picture [53][54][55][56] which describes the effect of the Gribov problem on correlation functions. Although originally analyzed in a semi-perturbative way, this can be generalized to a fully nonperturbative analysis [57][58][59]. Choosing a finite value of G(0), a family of so-called decoupling solutions is obtained. Their characteristic feature is an infrared (IR) finite gluon propagator [40,49,[60][61][62][63][64].
All other dressing functions are IR finite or logarithmically divergent. This can be shown analytically [49,51] but was also seen in various numerical calculations, e.g., [7, 30-32, 35, 41, 65-67]. The decoupling type of solutions can be accommodated in the Gribov-Zwanziger picture by taking into account certain condensates [61,68,69]. In addition, the finite value of the gluon propagator at zero momentum has motivated effective descriptions of correlation functions using massive extensions of Yang-Mills theory [70][71][72][73][74][75][76][77], e.g., based on the Curci-Ferrari model [70][71][72][73]78]. The origin of different solutions has been attributed to the Gribov problem which refers to the fact that the perturbative definition of the Landau gauge is only an incomplete gauge fixing [53,79,80]. Indeed, lattice calculations have seen variations of the correlation functions depending on the details of the employed gauge fixing algorithm [81][82][83][84][85][86][87][88][89]. However, a one-to-one correspondence between functional and lattice prescriptions has not been found yet, see, e.g. [89,90]. Here, the space of solutions is scanned by varying G(0). The scaling solution corresponds to the limit of an infinite gluon mass with G(0) → ∞. Lowering, G(0), also the gluon mass decreases. As discussed in detail in Ref. [41], there is a critical value which separates solutions into confined and Higgs-type classes. The former class is characterized by a maximum in the gluon propagator at non-vanishing momenta. For the Higgs-type solutions, the maximum is at zero. The existence of a maximum at nonvanishing momenta does not only lead to positivity violation of the propagator [15], it also reduces the spectral dimension, viz., the dimension felt by the propagator [91], from four to one. The calculations in this work span solutions from scaling to the boundary of the Higgs-branch. The gluon propagator is renormalized via momentum subtraction as well. However, there is an additional complication that haunted calculations of the gluon propagator for a long time, namely the appearance of quadratic divergences, see, e.g., [33][34][35][92][93][94][95][96][97]. Their origin lies in the breaking of gauge covariance by the regularization procedure. Methods to subtract these divergences are either limited in their practical applicability, for instance, because they require an exact knowledge of the vertices, or they introduce a new parameter. For an overview see Ref. [97]. The new parameter is a manifestation of the fact that the regularization scheme breaks gauge covariance. Consequently, a counter term for the gluon mass is no longer forbidden and can be used to renormalize the quadratic divergence [98]. However, the corresponding renormalization condition is not fixed and one would expect that results depend on what value one chooses. Indeed, calculations performed up to now using this renormalization scheme [13,96,99] saw such a behavior. The dependence on this parameter is also studied within the present truncation scheme in Sec. V E 2.
The gluon propagator renormalization is realized as follows [6,96]. We write the renormalized gluon propagator DSE as where Σ Z is the selfenergy, Z 3 the gluon wave function renormalization constant and C sub the renormalization constant for the quadratic divergences. One can determine the renormalization constant Z 3 in a standard MOM scheme by choosing a fixed value for Z(x s ): The term C sub is determined by demanding that the propagator D(x) = Z(x)/x has a fixed value at x m : In summary, we need to specify the renormalization conditions D(x m ) and Z(x s ) with their corresponding renormalization points x m and x s , respectively. In practical calculations, it turns out to be most convenient to choose x m at the lowest calculated point and x s in the perturbative regime. It should be noted that D(x m ) has a mass dimension of −2. The influence of choosing different val- For the scaling solution, the gluon propagator does not approach a fixed value at zero momentum but vanishes as D(x) = d gl x δ gl −1 . Thus, one cannot choose D(x m ) freely. On the other hand, it is known that on the righthand side of the gluon propagator DSE, the ghost loop is responsible for the IR suppression [33,42]. This does not change when two-loop diagrams are included [47][48][49][50][51] and leads to an analytic result for the IR fixed point of the MiniMOM coupling, given by by comparing the IR dominant diagrams of the ghost and gluon propagator DSEs [33,[42][43][44]. Via this relation, the value for D(x m ) can be determined from the ghost dressing function as explained in Appendix B.
The renormalization of the ghost-gluon vertex is trivial, as it is finite in Landau gauge [100]. Numerically, it turned out to be advantageous to increase the internally used cutoff for the vertex compared to the divergent quantities. This reduces subleading cutoff dependencies which cannot be absorbed in the renormalization process as for other divergent quantities. The employed renormalization scheme corresponds to the MiniMOM [33,101] or Taylor scheme [102] supplemented by the second renormalization condition for the gluon propagator. This scheme fixes all the renormalization constants via momentum subtraction for the propagators and a minimal subtraction for the ghost-gluon vertex, hence its renormalization constant is Z 1 = 1. The threeand four-gluon vertex renormalization constants, in turn, are fixed by Slavnov-Taylor identities (STIs). Since they are a consequence of the gauge invariance of the path integral, the STIs are modified by the presence of the ultraviolet (UV) cutoff in the numerical calculations. Thus, the original STIs relating the renormalization constants of the three-gluon and four-gluon vertices with those of the propagators and the ghost-gluon vertex, respectively, are only fulfilled approximately and they are not used for the renormalization of the respective vertices. Instead, a momentum dependent form of the STIs is used which contains the longitudinal pieces of the vertices, see also Ref. [41] for more details. Working in Landau gauge, the longitudinal pieces are not calculated here, but in the perturbative regime the longitudinal and transverse pieces agree. This entails that perturbatively the couplings of all vertices must agree. These running couplings are defined as follows [7,46,103]: The vertex dressing functions are taken at the symmetric points. 4 The perturbative equivalence of the couplings is then used to calculate a renormalization condition for the vertices from Note that this only ensures that the couplings agree at the renormalization point p r . Agreement over a wide range of momenta, as demanded by the STIs, is not guaranteed. This will be discussed further in Sec. V E 1.
The renormalization constants also appear in loop diagrams that feature bare vertices. In calculations using models for dressed vertices, they are often absorbed in the models which thus become cutoff dependent, see the discussion above. In the absence of such models as here, it is important to include the renormalization constants properly, as this is also important for the correct perturbative resummation of diagrams. It has to be noted that for the renormalization constants in the loop diagrams the values obtained from the STIs were used. In the iterative process, the renormalization constants need to be updated as well, and only once the final self-consistent solution is obtained, the renormalization constants have obtained values that balance all equations. This introduces a certain destructive element and it turned out to be most stable to use the renormalization constants from the STIs in these cases. For the final solution, the values calculated from the two methods agree better than 1 %. In addition, the renormalization constants of the propagators also agree with the inverse dressing functions at the cutoff to the same degree. In previous calculations, these quantities typically deviated by a few percent.

V. RESULTS
This section contains the results for the propagators and vertices. Furthermore, self-tests of the results are discussed and performed. All calculations were done for the gauge group SU (3). For reference, plots contain also lattice data where available. A more detailed comparison to lattice and FRG data is done in Sec. VI.

A. Propagators
As discussed in Sec. IV, the equations allow a continuum of solutions by choosing different values for the ghost dressing function in the IR. For the gluon and ghost propagators, results are shown in Figs. 8 and 9, respectively. The spectrum of solutions is represented by the gray band which is bounded by the scaling solution (dashed line). The other bound is a decoupling solution close to the Higgs-branch of solutions (black continuous line), viz. the solution with the maximum at the smallest non-vanishing momentum found in the set of calculated decoupling solutions. As a consequence, the maximum is extremely shallow and can only be seen directly in the data and not in the plots. Intermediate decoupling solutions are shown as thin red lines. All solutions agree roughly down to 2 GeV. Below, the ghost dressing function, shown in Fig. 9, bends up and becomes finite for all solutions except the scaling solu-  [43,44]. Extracting momentum dependent IR exponents, viz. calculating the power laws for small momentum intervals, one observes that the IR power laws of the scaling solution are reached also for the decoupling solutions close to the scaling solution before they turn to their final values δ gh = 0 and δ gl = 1. The comparison with lattice results does not show a candidate solution that agrees well in all quantities. The best overall agreement is obtained with the solution close to the Higgs-type branch. However, this agreement is not perfect and one cannot establish a one-to-one correspondence between lattice and functional results. This was also observed for the corresponding FRG solutions [41]. Both ghost and gluon propagators show, though, that the lattice results are closer to the Higgs-branch than to the scaling solution. This also explains why the maximum in the gluon propagator is difficult to observe in lattice results for four dimensions [62,104,105]. In general, however, there are additional arguments for the existence of a maximum from both continuum analyses, e.g., [91,106] and from three-dimensional lattice calculations, e.g., [88,107,108]. The observed maximum in the gluon propagator directly leads to a violation of positivity of the spectral function, see, e.g., [86,91,109,110]. This is reflected in the Schwinger function calculated from the propagator as In Fig. 10, the Schwinger function for different solutions is shown. For the scaling solution, the violation of positivity around 1 fm is most pronounced, but also the other solutions become clearly negative. It should be noted, though, that in the Landau gauge positivity is violated already perturbatively [15,111].

B. Ghost-gluon vertex
Results for the ghost-gluon vertex are shown in Fig. 11. Due to the comparatively large angular dependence of this vertex, it is worthwhile to compare different kinematic configurations. Again, a selection of the family of solutions is shown. At low momenta, the vertex dressing function becomes constant. However, while it becomes one for all decoupling solutions, it settles at a value larger than one for the scaling solution. This difference can be understood analytically from the IR behavior of the propagators and vertices in the integrals. It is known that with the boundary condition of a divergent ghost dressing function [43,44], vertex dressing functions behave like (p 2 ) (−m+n/2)κ in the IR [46], where m/n is the number of gluon/ghost legs. This behavior is induced by diagrams with a bare ghost-gluon vertex in the DSE [47,50] and, consequently, there is an IR constant contribution from these diagrams that change the total IR value away from one given by the tree-level. For a decoupling solution, on the other hand, one can show that all diagrams are IR suppressed [49] and thus the tree-level value is not altered. Numerically, this behavior was already confirmed for scaling [112] 5 and decoupling [35,113]. It should be noted that for the scaling solution the IR limit of the ghost-gluon vertex depends on the direction from which one approaches zero momentum. When comparing the ghost-gluon vertex results to lattice results, it becomes evident that there is a mismatch in the scale. Most likely, this is not due to comparing to SU (2) lattice data, since the differences between SU (2) and SU (3) are only weak [114,115]. However, the position of the bump is found at lower momenta in all continuum results [35,65,67,116,117]. As can be seen, the scaling solution is the solution with the bump at lowest momenta. Moving away from the scaling solution, the position of the bump moves to higher momenta. In the bottom right plot, the full angular dependence for the lowest decoupling solution, corresponding to the black continuous lines in the other plots, is shown together with lattice data for various kinematics. For a meaningful comparison, the data is plotted over the variable S 0 , see Eq. (15). For this single solution, the agreement with the lattice data is better than for the other solutions, but there is room for improvement. From all solutions, though, this is the one closest to minimal Landau gauge as employed on the lattice. This hypothesis is supported by similar observations for the propagators, see Sec. V A.

C. Three-gluon vertex
The three-gluon vertex shows a remarkable small angular dependence as illustrated previously in Fig. 5 where the dependence on the momentum scale S 0 is shown and the small band represents the variation of the dressing function with the variables a and s. As a consequence, it is sufficient to compare with other results for one momentum configuration only. For the sake of comparison, all results were renormalized to one at 5 GeV. Fig. 12 shows the family of solutions and a comparison to var- ious lattice results for the symmetric configuration. All solutions cross zero and diverge; logarithmically for decoupling solutions and like (p 2 ) −3κ for the scaling solution. The position of the zero crossing correlates with the position of the maximum of the gluon propagator as argued for in Ref. [106]. Extrapolating the position of the zero crossing for the case of a gluon propagator with no maximum is compatible with the vanishing of the zero crossing. Again we find that the decoupling solution close to the Higgs branch agrees best with lattice results. It should be noted that the three-gluon vertex has a non-trivial IR structure in the scaling case [120], which, however, is not visible here due to the employed projection onto the tree-level tensor.

D. Four-gluon vertex
The four-gluon vertex is calculated with the singletdoublet approximation for the kinematics and only one dressing function as discussed in Sec. II. The result is shown in Fig. 13 where one can see that the angular dependence is stronger than for the three-gluon vertex. However, the breadth of the dressing function around 10 GeV comes from a few outlying points whereas the majority clusters in a much smaller band. This is indicated by plotting besides the band also individual points in Fig. 13. The qualitative behavior is in agreement with previous calculations [32,41,66,121,122]. In the IR, the vertex diverges like (p 2 ) −4κ for the scaling solution. For the decoupling solution, no logarithmic IR divergence is found for the tree-level dressing function in agreement with previous calculations [32,66]. It should be noted that the individual ghost box diagrams are logarithmically divergent, but the symmetrization of the equation cancels this effect. A small number of other tensors has already been explored [32,66]. Interestingly, in that case logarithmic divergences were found [32,66], but the cor-responding dressing functions are smaller than the treelevel dressing function. There are 5 × 41 tensors in total, and the question of their quantitative importance remains open. A remarkable feature of the four-gluon vertex is that the zero momentum limit depends on the kinematic configuration for which it is approached, as can be seen by the finite width of the band in the IR in the left plot of Fig. 13. This is similar to the ghost-gluon vertex for the scaling solution. For the three-gluon vertex, any such dependence is at best too small to allow any conclusions. Such an IR irregularity is interesting insofar, as it was argued in Ref. [41] that the creation of the gluon mass gap for decoupling solutions requires vertices that show an irregular IR behavior. One should be careful, though, before drawing any final conclusions, because currently the four-gluon vertex is calculated with the most severe approximations of all quantities.

E. Self-tests
Assessing the reliability of results from functional equations is one of the main challenges of this approach. Beyond some trivial self-consistency checks, as, for example, the correct asymptotic IR and UV behaviors, a standard test is to compare to other results, typically from lattice calculations. However, this has two drawbacks: First, one relies on the availability of other results. Second, even more importantly, such comparisons should not necessarily be interpreted too literally on a quantitative level. For instance, the renormalization schemes are different for lattice and functional methods. Actually, also within the latter category care has to be taken when comparing results from different calculations. While for DSE or nPI calculations typically a MOM scheme is employed, the renormalization scheme is implicit for FRG calculations. Furthermore, the correspondence of mem- bers of the family of solutions between different methods is not clear due to the different ways they are realized for different methods, e.g., by changing the ghost dressing function at zero momentum, by choosing a UV mass parameter or by a selection algorithm of Gribov copies. In this situation, any way of checking the results without requiring external input is welcome. Here, two possibilities are explored. The first relies on the different vertex couplings and the requirement that they should agree perturbatively. The second self-test checks the irrelevance of an unphysical parameter. In addition, the application of the present results to the calculation of a gauge invariant object is discussed.

Couplings
As discussed in Sec. IV, the couplings extracted from the different vertices should agree in the perturbative regime. This turns out to be a nontrivial requirement and can only be achieved when dynamically including all required quantities. Such a computation was performed in Ref. [41] and the couplings were found to agree down to a few GeV. There, one can also find a comparison between different previous calculations which calculated different subsets of the system considered here. The comparison reveals that such isolated calculations are insufficient as deviations between the couplings were found. In this work, for the first time, all these quantities are calculated from their equations of motion jointly and selfconsistently which leads to a good agreement between the various couplings. The couplings from the ghost-gluon, three-gluon and four-gluon vertices, defined in Eq. (52), are shown in Fig. 14. They show good agreement down to 3 GeV. The four-gluon vertex coupling deviates earlier from the ghost-gluon vertex coupling than that of the three-gluon vertex. In Sec. VI, the couplings are also compared to results from the FRG. It has to be stressed that the employed renormalization conditions for the vertex, see Sec. IV, are not sufficient to obtain this degree of agreement, because they only force the couplings to agree at one specific point and do not fix their running. Thus, the good agreement over several orders of magnitude comes from the dynamics of the equations.
The couplings have two convenient properties: They are independent of the renormalization point and their runnings relate directly to the scale of Yang-Mills theory Λ YM . The former property makes the couplings a useful quantity for comparisons between different calculations. The latter property could be used to carry over the scale from lattice simulations. However, in the perturbative regime, for which dedicated lattice simulations exist to calculate the scale parameter of Yang-Mills theory Λ YM [102,123], the slow logarithmic running of the coupling does not allow a precise determination of the scale for our purposes. Thus, the scale fixing is done in the nonperturbative regime where it is most convenient to fix the scale via the position of the bump of the gluon dressing function for which p 0 = 0.97 GeV was chosen. Even though the different solutions differ already in this regime, this procedure is accurate enough as can be seen by the good agreement of the different solutions above 1 GeV.
Having fixed the scale like this, another nontrivial check is to compare the couplings in the perturbative regime. For this, the MiniMOM coupling defined in Eq. (50) is most convenient as we can compare directly to high momentum results of Ref. [123]. As reference points, we take the values p 2 = 178 GeV 2 and 1785 GeV 2 , which correspond to 1, 000 and 10, 000 in units of the Sommer scale r 0 . The results for the coupling from the lattice calculation are roughly 0.142 and 0.107, whereas the DSE cal- culation yields 0.144 and 0.107, respectively. This good agreement further supports the quantitative reliability of the results and it should be stressed that there is no parameter to tune these values. The use of a hard UV cutoff entails that gauge covariance is broken and the STIs need to be modified. As a consequence, the momentum independent STIs given in Eq. (51) are not fulfilled exactly. However, observing the good agreement between different couplings is a manifestation of the restoration of gauge covariance. This can also be tested by comparing the values for the vertex renormalization constants calculated in their renormalization as discussed in Sec. IV with the values obtained via the STIs given in Eq. (51) from the propagator renormalization constants. The differences are approximately 1 % for Z 1 and 0.2 % for Z 4 which is an improvement over previous calculations where these differences were at the order of a few percent.

Quadratic divergences
Another consequence of the hard UV cutoff is the appearance of quadratic divergences, see the discussion in Sec. IV. To deal with them, a second renormalization condition D(x m ) was introduced which is given by the value of gluon propagator at a given IR scale x m . One could expect that this is similar to how the problem of quadratic divergences is handled in the flow equation of the gluon propagator where the value of the bare gluon mass is fine-tuned and leads to different solutions in the IR [41]. Indeed, a first test for three-dimensional Yang-Mills theory using a simple truncation showed that the gluon mass renormalization can be used to obtain a family of solutions [99]. However, the more sophisticated truncation of the present work sheds new light on this method to renormalize the gluon propagator DSE. As it turns out, the results are basically independent of the new renormalization parameter D(x m ). Necessarily, since this parameter is a dimensionful quantity, the scale setting factor changes, but the solutions are equal, as can be seen by comparing any of the running couplings, which are RG invariant. The propagators and vertices themselves are merely rescaled as required by multiplicative renormalizability. Taking this finite renormalization into account, quantities can be compared for different values of D(x m ). For illustration, the gluon propagator and its dressing function are shown in Fig. 15 for the three values D(x m ) = 10, 12 and 14 (in internal units). The plots show the data in physical units renormalized at 6 GeV. The small differences below 2 GeV could be both numeric or truncation artifacts. Since quadratic divergences are a manifestation of broken gauge invariance, the fact that they can be unambiguously subtracted seems naturally related to the good agreement of the different couplings which is a consequence of gauge invariance itself. In turn, the small remaining differences might reflect the small deviations of the couplings discussed in Sec. V E 1. However, quantitatively they are so small that they do not matter for physical applications like the calculation of glueballs discussed in Sec. V E 3. The differences in the other dressing functions are even smaller than for the gluon propagator.
The absence of a parameter related to the removal of quadratic divergences is a novel feature encountered for the first time within a DSE calculation. The original idea that both the FRG and DSEs can create a family of solutions by varying a mass parameter related to the removal of quadratic divergences might have seemed appealing [99], but the fact that these subtraction methods are not equivalent settles the problem of having a second parameter for DSEs, the ghost renormalization condition, that fulfills the same role. Thus, we are finally in the position to subtract quadratic divergences unambiguously which is required to provide quantitatively reliable results.

Glueballs
An ultimate test of the quantitative correctness of the gauge-dependent propagators and vertices obtained here is employing them for the calculation of gauge-invariant quantities. The obvious candidates for Yang-Mills theory are glueballs, i.e., bound states of gluons. The simplest case are scalar and pseudoscalar glueballs which can be calculated from Bethe-Salpeter equations as twogluon bound states [5,[124][125][126]. In Ref. [5], the kernels of these Bethe-Salpeter equations were constructed from the same truncation of the 3PI effective action as here. This provides a fully self-consistent setup. The masses of the ground state and two excited states were then calculated using the propagators and vertices obtained here as input [5]. For the scalar glueball, the ground state mass was calculated as 1.75 ± 0.12 GeV and for the pseudoscalar glueball as 2.44 ± 0.17 GeV. These results agree well with corresponding lattice results of 1.73 ± 0.13 GeV and 2.59 ± 0.17 GeV, respectively [127]. Also the first excited states agree well: The bound state calculation yields 2.43 ± 0.2 GeV and 3.65 ± 0.11 GeV for the scalar and pseudoscalar masses, respectively, while on the lattice 2.67 ± 0.31 GeV and 3.64 ± 0.24 GeV are found. Second excited states in each channel were also calculated. The calculation of gauge-invariant quantities conveniently also provides the means to test if different solutions from the family of solutions are indeed physically equivalent. To this end, the calculation of the glueball masses was repeated with different sets of propagators and vertices. As no difference in the bound state masses was found, this provides another piece of evidence that different solutions just represent different nonperturbative completions of the perturbative Landau gauge. Note that if the glueball masses varied with the input, one could not decide whether the employed truncation is insufficient or if the solutions really differ physically.
To test the impact of the self-consistency of the input, also mixed solutions were used, viz., propagators and vertices were taken from different solutions. It should be noted that mixing solutions like this requires to rescale them properly, i.e., fix them to the same scale, as otherwise the different scales of the propagators and vertices destroy the consistency automatically. For such mixed input, the bound state masses maintain the correct hierarchy but deviate when the input solutions are taken from members of the family of solutions which are close to each other. If the mixed solutions are too far apart, though, the hierarchy of the masses gets lost and no sensible solutions can be obtained anymore.
As a final test, the solution from the system with the ghost-gluon vertex DSE instead of its EOM, discussed in App. A, was used. The glueball masses were equal within the given errors.

VI. COMPARISON WITH OTHER RESULTS
The results for all correlation functions are compared in this section to results from the lattice and the FRG. In general, good agreement is found. The purpose here is to identify and highlight the remaining differences. The system of flow equations is formally truncated in the same way as here by including all primitively divergent correlation functions and neglecting the others. Due to the different structure of the equations, the truncations are not equal, though. For example, the two-ghost-two-gluon and the four-ghost vertices appear in the propagator flow equations. The corresponding diagrams were dropped in Ref. [41] but investigated for three dimensions [128]. The FRG results shown here correspond to the '1D momentum dependent' approximation of Ref. [41], i.e., the three-point functions are calculated only at the symmetric point. As is clear from the small angular dependence of the three-gluon vertex, this is a good approximation in this case. For the ghost-gluon vertex it is more severe, but the impact is still small [41]. The four-gluon vertex is solved for the FRG for two kinematic configurations.
In the plots the DSE and FRG results are distinguished by colors. The scaling solution is represented by a dashed line, the decoupling solutions by continuous ones. For the comparison, all results were renormalized such that they agree in the region around 5 GeV.
The propagator results are compared in Figs. 16 and 17. From the ghost dressing function and the gluon propagator one can see that for both functional methods the agreement with lattice data is best for a decoupling solution with a small gluon mass gap. It should be noted that the gluon dressing function is slightly affected by the '1D' approximation of the FRG calculation, but for the purpose of the comparison here this effect is small enough. In general, it is not possible to find a solution for which both propagator dressing functions agree between the FRG and DSE results. The cleanest compar-ison can be done for the scaling solution, because there is only one instance. However, it can not be expected that there is a clean one-to-one correspondence of single members of the family of solutions, because the similarities of the truncations are only superficial. For example, the role of the tadpole diagram in the FRG calculation is partially played by the one-loop gluon, the sunset and the squint diagrams in the DSE. Thus, the treatment of the four-gluon vertex affects both systems differently. For the ghost-gluon vertex, the comparison is shown in Fig. 18. Both functional results are very similar, and one should keep in mind that the '1D' approximation is most severe for this quantity. In particular, the FRG results show a similar position of the maximum as the DSE results which differ from the lattice results. As mentioned above, this is a generally observed difference between lattice and continuum results. In addition, the shift of the position of the maximum from lower to higher momenta from the scaling solution through the decoupling solutions is similar for both functional methods. For the three-gluon vertex, shown in Fig. 19, one finds an interesting difference between the DSE and FRG results. The data is renormalized at 5 GeV. As can be seen, the results deviate in the region below 2 GeV. The FRG results have larger dressing functions at lower momenta. If one tries to match the results to the lattice data by adjusting the normalization, the DSE results show a larger deviation than the FRG results. These differences are most likely due to the different truncations. As was discussed in Sec. III C, there are also differences between a two-loop DSE and three-loop truncated 3PI EOM. They are small, but the two-loop DSE results are also larger at low momenta. Also for the four-gluon vertex deviations are observed, especially below 3 GeV. It can be seen in the present (see Fig. 13) and previous DSE calculations [32] that there is a large angular dependence compared to the three-gluon vertex. Not taking this into account most likely affects  [41] and SU (2) lattice data [90] (details as explained in Fig. 11).  [41] and lattice data [118,119] (details as explained in Fig. 12).
the FRG calculation and also here only three out of six kinematic variables were taken into account. This may explain at least a part of the quantitative differences. As final quantities, the vertex couplings given in Eq. (52) are compared in Fig. 21. The running of the couplings agree quite well, but it is found that they deviate by an overall factor. This is taken into account by rescaling the FRG couplings such that they agree with the DSE results at 10 GeV. For reference, the ghost-gluon vertex coupling is also shown without rescaling. Whether this difference is entirely due to differences in the renormalization schemes or due to other sources still has to be investigated.
The largest difference is seen in the four-gluon vertex coupling. This is most likely a direct effect of the different approximations. It should be noted that the fourgluon vertex is here calculated from its DSE and not an EOM, since at least the 4PI effective action is re- quired for an EOM. It is also confirmed that the effect of kinematic approximations cannot be solely responsible for this deviation, as a calculation with the vertex restricted to the symmetric point yields a very similar coupling. The ghost-gluon and three-gluon vertex couplings agree very well down to 4 GeV, for the FRG even further to approximately 2.5 GeV.
In general, the agreement between the vertex couplings is a desirable feature of any calculation of correlation functions which can only be realized in self-consistent truncations. In many previous studies, this is not found and the runnings of the couplings disagree already in the perturbative regime. To what extent the agreement needs to be realized is not fully clear yet and might depend on the specific problem. For example, for the calculation of glueball masses the present degree of agreement was found to be sufficient [5].

VII. DISCUSSION
The present truncation includes all primitively divergent correlation functions of Yang-Mills theory. While the propagators and the ghost-gluon vertex are calculated entirely, for the three-and four-gluon vertices approximations of their full forms were made. The approximation for the three-gluon vertex consists of taking into account only one out of four transverse dressing functions. It was found previously that the basis can be chosen such that the dressing function of the tree-level tensor dominates [7]. Indirectly, the results here corroborate that, but it remains to be seen what effect additional, albeit small contributions of the three-gluon vertex can have.
In particular, the gluon propagator DSE is quite sensible to changes in the midmomentum regime, simply because the gluon propagator is the inverse of a sum of diagrams. Small changes around the peak of the gluon dressing function can thus have a large effect, which is one of the reasons why the subtraction of quadratic divergences needs to be done with high enough precision.
For the four-gluon vertex the same approximation was made with regard to the tensor basis. In this case, no calculation using the full tensor basis exists. Only projections onto some tensors were tested and found to be smaller than the dressing function of the tree-level tensor [32,66]. The four-gluon vertex enters via the sunset diagram in the gluon propagator DSE. This softens its possible effect on the gluon propagator, because this diagram is subleading compared to the one-loop gluon diagram and the squint diagram through which the three-gluon vertex enters [29]. However, as discussed in Sec. III, the presence of the dressed four-gluon vertex has within the present truncation a stabilizing effect. A second approximation for the four-gluon vertex is made with respect to its kinematic dependence. Instead of six kinematic variables, only three are taken into account to reduce the computational costs. A qualitative comparison with the results of Ref. [32], where the full kinematic dependence was considered, shows a similarly large angular dependence. For a detailed estimation of the induced error, though, a dedicated calculation with all six kinematic variables is required.
To assess the truncation error, several tests of the results were performed. One crucial test is the good agreement of the vertex coupling down to a few GeV. This was also found in FRG calculations [41], for which the agreement is even better than observed here. The reason is most likely that here the vertices are not treated on the same footing: The three-point functions are calculated from their EOMs while for the four-gluon vertex its DSE is used. Thus, the agreement between the ghost-gluon and three-gluon vertex couplings is better than that of the four-gluon vertex coupling with either of them. However, the difference is small enough to not affect the calculation of glueballs and the handling of quadratic divergences.
It is an open question in what cases the found agreement is sufficient and when an improvement will be required.
A second test concerns the problem of quadratic divergences. The treatment of these divergences in previous calculations of the gluon propagator DSE should be considered part of the truncation. Here, on the other hand, the solution was shown to be independent of the parameter introduced by the subtraction of the divergences. Strictly speaking, it cannot be excluded that the disappearance of this dependence is a lucky coincidence. However, together with the other findings, it seems likely that the problem is solved. Closely related to this is the good agreement of the vertex couplings mentioned above, as both issues are related to the breaking of gauge invariance. A successful solution to both problems hints at an effective restoration of it. However, in both cases there is still room for improvement as illustrated by the remaining tiny dependence on the subtraction parameter, illustrated in Fig. 15, and the not yet perfect agreement of the vertex couplings down to scales of 2 − 3 GeV. It has to be stressed that, given the good quantitative agreement with lattice results, the errors have reached a new qualitative level compared to previous results from equations of motion. As a consequence, the application of these results to other problems is promising. The calculation of glueballs discussed in Sec. V E 3 is one example where this was successfully realized. An important issue for all functional calculations is the stability of an employed truncation. Naturally, one cannot exclude the possibility that two or more individual effects somehow cancel each other. Hence, adding only one of them would lead to wrong conclusions. An example of such an effect was observed in three-dimensional Yang-Mills theory where it was found that the effect of including a certain correlation function into a system of flow equations was counteracted by the inclusion of another one [128]. For the truncation employed here, the situation with regard to stability of the truncation can be summarized as follows. The impact of neglected higher correlation functions was partially tested already in the past. This concerns diagrams with the two-ghost-two-gluon and four-ghost vertices [13], while correlation functions beyond four-point functions are currently untested. For the three-gluon vertex, it was shown that a one-loop truncation of its DSE leads to a qualitative defect and thus this truncation can be considered insufficient. On the other hand, the EOM from the three-loop truncated 3PI effective action and a DSE including most two-loop diagrams yield very similar results, which are in addition close to lattice and FRG results. This indicates that the three-gluon vertex is described quite well, but some small differences remain which could be corrected by extensions of the present setup. The natural candidate would be the EOM from a 4PI effective action to include a dressed four-gluon vertex. For the ghost-gluon vertex the situation is different. Using EOM or the DSE yields different results, but at the same time the other correlation functions remain basically unaffected, see Appendix A. This could reflect a general stability of the system of equations.
It should be stressed that this is by no means a trivial observation, as correlation functions not only react quantitatively to changes in other ones but often the convergence of the system itself is jeopardized by shortcomings in the truncation. Two concrete examples are the gluon propagator and the three-gluon vertex. It is found that the dressing function of the former cannot have a much larger maximum, as then the DSE does not converge anymore. In a sense, the gluon propagator solution exists close the border of possible solutions. For the three-gluon vertex, a balance between the gluon triangle and the other gluonic diagrams is required, as a too strong gluon triangle, as caused by a too large gluon propagator, prevents the equation from converging.

VIII. CONCLUSIONS
The parameter-free calculation of the correlation functions of QCD is a necessary step for functional equations to make the transition from providing effective descriptions based on models fixed by phenomenology to establishing them as a first principles method. Solving a given set of equations is one challenge in this regard.
Proving that the results are indeed sufficiently close to the correct solution is a second. In this work, all primitively divergent correlation functions of Yang-Mills theory were solved for the first time simultaneously from a coupled system of equations of motion. To address the second challenge, several tests were performed that show the qualitative and quantitative improvements compared to previous calculations. In general, good quantitative agreement with results from other methods, namely lattice simulations and the FRG, is found which provides additional reassurance about the reliability of the results. At the same time, with this level of precision one can now focus on the remaining differences. Albeit being small, they might play an important role in the further development of truncations. In partic-ular, one should learn which deviations are relevant and which are not.
Several extensions of the this work present itself. Naturally, it would be interesting to bring calculations with other gauges to the same level of truncation. In fact, in some cases, like the maximally Abelian gauge, this seems like a requirement for numeric calculations anyway, because no simpler working truncation was found yet, see Refs. [6,50,129] for details. In other cases, though, solutions of simple truncations exist. An example are linear covariant gauges [130,131], for which, however, the treatment of the longitudinal sector will be an additional challenge. Another example is the Coulomb gauge, for which the Hamiltonian approach provides an alternative functional approach that was already successfully employed in the past, e.g., [132][133][134][135][136][137][138].
For Yang-Mills theory in the Landau gauge, it would be interesting to consider the 4PI effective action and its equations of motion to put all correlation functions on the same footing. This should then further improve the agreement of the four-gluon vertex coupling with the other couplings. Physically, the interesting extension is the inclusion of quarks. The quark sector was already calculated with the 3PI effective action, but the Yang-Mills propagators were used as fixed input [139]. From the delicate balancing of the equations in the Yang-Mills sector found here and from FRG calculations [140], it can by no means be expected to be trivial to include quarks. Nevertheless, such calculations will be necessary to describe full QCD in a self-contained way with this method.
A successful inclusion of fermions would also pave the way for addressing additional interesting questions with functional methods in a qualitatively new way. For example, technicolor scenarios could be probed or the effects of temperature and density in QCD with a largely reduced or even eliminated need for modeling. Such questions are difficult to answer with any method, but the recent progress with functional methods leads to new perspectives. As discussed in Sec. III C, the EOM from the three-loop truncated 3PI effective action is used for the ghost-gluon vertex. Since its DSE looks very similar, one might wonder if this equation provides an equally good description. To this end, the same system of equations was solved but with the ghost-gluon vertex EOM replaced by its DSE. A comparison between the DSE and EOM solutions is shown in Fig. 22 where it can be seen that the strength of the vertex dressing function changes roughly by 10 %. From the lower right plot one might even think that the DSE solution agrees better with the lattice results than the EOM solution. However, as can be seen from specific kinematic configurations, as shown in the other panels, this is not the case. An interesting finding is that the dressing functions of the other correction functions, which are shown in Fig.  23, do not change despite the differences in the ghostgluon vertex. This independence of all quantities from the ghost-gluon vertex is rather unexpected. For exam-ple, the ghost propagator is known to react to the used ghost-gluon vertex input, e.g., [35,65,142]. The absence of such a dependence could be another indication of the stability of the presented solution for the full system. and gluon DSEs can be written as follows: I gh and I gh are functions of δ gh and δ gl which can be calculated analytically for the IR limit. From the two equations, one can infer that I gh = I gl which allows to calculate values for δ gh and δ gl . We use this equality in Eq. (B1b), to calculate d −1 gl . It remains to rewrite the equation in terms of the ghost self-energy Σ G (x) = g 2 N c d gh d gl x δ gh +δ gl I gh and the dressing functions: . (B3) Note that d −1 gl also appears on the right-hand side. The equation is trivially true for the solution of the system of equations, but not during the iteration process. Hence, to update d gl , Eq. (B3) is used and D(x m ) is calculated from it. Technically, it might be worthwhile to point out that in this case starting guesses for the dressing functions were used for the iteration that were obtained using the IR fixed point value and a very small relaxation parameter for the gluon. This led to a stable solution and no further study of the stability of this procedure for different starting guesses was performed.

Appendix C: Numerical implementation
This section describes some numerical details of the calculations. As basis, CrasyDSE [27] was used that provides basic C++ routines for integration and interpolation.
• The three-dimensional interpolation for the vertices uses cubic spline interpolation for S 0 , linear interpolation for ρ and trigonometric interpolation for η.
• Starting conditions: For most calculations, results from previous calculations were used, as this turned out to be more convenient than using generic starting ansaetze.
• Iterative process: The calculation used three levels of iteration. In the innermost loop, equations were iterated by themselves. However, typically, only for the ghost propagator more then one iteration step was done. All equations were then iterated consecutively (meta iteration), starting with the propagators. After at most 10 iterations, the renormalization constants were updated and the meta iteration started again (super iteration). Normally, only the ghost propagator equation required relaxation.
• Integration: Standard Gauss-Legendre integration was used for all integrals. The integration intervals for the propagators were split at the external points, as explained, e.g., in Ref. [143]. For the vertices, also appropriate splittings for the intervals were performed.
• Extrapolation: The propagators were extrapolated in the IR and UV by the known analytic forms. For the ghost-gluon vertex, the boundary values were employed. The three-and four-gluon vertices were extrapolated in the IR by their boundary values. It was tested that this does not affect the results. In the UV, STI motivated extrapolation functions proportional to Z(x)/G(x) and Z(x)/G(x) 2 , respectively, were used.
• Scale setting: The coupling in the calculation was set to α(µ 2 ) = 0.05. This determines the scale, which is determined a posteriori. However, instead of determining the physical value of µ from the coupling α(µ 2 ) = 0.05, the scale is fixed by putting the maximum of the gluon dressing function at p 0 = 0.97 GeV.