On mass generation in Landau-gauge Yang-Mills theory

A longstanding question in QCD is the origin of the mass gap in the Yang-Mills sector of QCD, i.e., QCD without quarks. In Landau gauge QCD this mass gap, and hence confinement, is encoded in a mass gap of the gluon propagator, which is found both in lattice simulations and with functional approaches. While functional methods are well suited to unravel the mechanism behind the generation of the mass gap, a fully satisfactory answer has not yet been found. In this work we solve the coupled Dyson-Schwinger equations for the ghost propagator, gluon propagator and three-gluon vertex. We corroborate the findings of earlier works, namely that the mass gap generation is tied to the longitudinal projection of the gluon self-energy, which acts as an effective mass term in the equations. Because an explicit mass term is in conflict with gauge invariance, this leaves two possible scenarios: If it is viewed as an artifact, only the scaling solution survives; if it is dynamical, gauge invariance can only be preserved if there are longitudinal massless poles in either of the vertices. We find that there is indeed a massless pole in the ghost-gluon vertex, however in our approximation with the assumption of complete infrared dominance of the ghost this pole is only present for the scaling solution. We also put forward a possible mechanism that may reconcile the scaling solution, with an infrared dominance of the ghost, with the decoupling solutions based on longitudinal poles in the three-gluon vertex as seen in the PT-BFM scheme.


I. INTRODUCTION
One of the central open questions in strong interaction studies is the origin of mass generation in Quantum Chromodynamics (QCD). One aspect of the problem is the fact that the majority of the masses of light hadrons, and therefore the visible mass in the universe, must be generated in QCD because light current quarks only carry a small fraction of the mass of the proton. Mass generation in the quark sector is relatively well understood by now in terms of dynamical chiral symmetry breaking and the corresponding dynamical generation of a large quark mass at low momenta, which is seen in various nonperturbative approaches including lattice QCD [1][2][3][4] and functional methods such as Dyson-Schwinger equations (DSEs) [5,6] and the functional renormalization group (fRG) [7,8].
Another and perhaps more fundamental aspect is the emergence of a mass gap in pure Yang-Mills theory, i.e., QCD without quarks. This is tied to the open question of confinement, and the corresponding problem of mass generation in the Yang-Mills sector of QCD is much less understood.
In principle, the origin of mass generation is encoded in QCD's elementary n-point correlation functions. For Yang-Mills theory, these are the two-point functions (the gluon and ghost propagators), three-point functions (the three-gluon vertex and ghost-gluon vertex), four-point functions (e.g., the four-gluon vertex) and higher n-point functions. In particular, it is well-established by now that the massless pole in the perturbative gluon propagator cannot survive in nonperturbative calculations in general covariant gauges including Landau gauge. Moreover, the respective mass gap is directly related to confinement as shown in [9,10].
Possible mechanisms for the generation of the mass gap include the Kugo-Ojima confinement scenario [11], where the n-point functions scale with infrared (IR) power laws [12][13][14][15][16][17][18][19][20] as given in Eq. (4) below, a Schwinger mechanism for longitudinal correlation functions [21][22][23][24][25][26][27][28], and the related gluon condensation mechanism [29,30]. In all these scenarios, irregularities in longitudinal and/or transverse projections of the Yang-Mills vertices are triggered and required for the mass gap to be present. While the IR scaling of the correlation functions in the Kugo-Ojima scenario directly induces such irregularities by assuming the existence of BRST charges, these must be triggered explicitly in the Schwinger mechanism and gluon condensation scenarios.
The gluon propagator is parametrized as where are the transverse and longitudinal projection operators with respect to the four-momentum Q µ , ξ is the gauge parameter and ξ = 0 corresponds to Landau gauge. In linear covariant gauges, the longitudinal dressing function L(Q 2 ) = 1 is trivial due to gauge invariance, but for later purposes we keep it general in what follows. A basic mass parameter m 0 is defined by the value of the transverse part D(Q 2 ) = Z(Q 2 )/Q 2 at Q 2 → 0, The IR solutions seen in lattice QCD simulations do not show scaling as required in the Kugo-Ojima confinement scenario and have been called decoupling or massive solutions. In this case, the gluon propagator saturates at low momenta and becomes constant in the IR [31][32][33][34][35] like in Eq. (3). This solution is also obtained with DSE and fRG calculations, see e.g. [23,25,[36][37][38][39][40][41][42][43]. In addition, functional studies of correlation functions within approximations have revealed the existence of a family of these decoupling solutions, where the maximal decoupling solution is close to that found on the lattice. There are also indications for different decoupling solutions on the lattice depending on the IR details of the gauge-fixing procedure, which involves the removal of Gribov copies [44]. In the continuum limit this is a numerically very challenging problem which has not been overcome yet. It has been speculated that the emergence of a family of solutions may be due to an additional gauge fixing parameter in Landau gauge, see e.g. [45,46].
Correspondingly, for the decoupling solutions the transverse gluon dressing function Z(Q 2 ) vanishes like Z(Q 2 ) ∝ Q 2 and 1/Z(Q 2 ) ∝ m 2 0 /Q 2 has a singularity at the origin in Q 2 . This is shown in Fig. 1 and clearly differs from a QED-like behavior where the photon dressing function and not the propagator saturates at IR momenta. The origin, details and consequences of this feature are however not yet fully understood. A possible explanation that has been studied in the PT-BFM (Pinch Technique/Background-Field Method) framework is due to the Schwinger mechanism, which could generate longitudinally coupled massless poles in Yang-Mills vertices such as the three-gluon vertex and thereby induce such a behavior in the transverse part of the gluon propagator [21,22,[24][25][26][27][28].
The other endpoint of the family of decoupling solutions for m 0 → ∞ is the scaling solution [12][13][14][15][16][17][18][19][20] with the IR scaling Z(Q 2 → 0) ∝ (Q 2 ) 2κ , G(Q 2 → 0) ∝ (Q 2 ) −κ . (4) Here, G(Q 2 ) is the ghost dressing function and the IR exponent κ is typically of the order κ ∼ 0.6 depending on the truncation of the system. Such a behavior is consistent with the Kugo-Ojima confinement scenario based on the assumption of global BRST symmetry (existence of BRST charges); see [39,47] for detailed discussions. By contrast, the decoupling scenario entails In any case, the gluon propagator is neither gauge invariant nor renormalization-group invariant and thus its value (3) at vanishing momentum hardly defines a mass gap without further specifications. Instead, a gluon mass gap is best defined as the (spatial) transverse correlation length of this correlation function through the screening mass, see e.g. [48] for the finite temperature version. It can be extracted from the gluon propagator D µν (Q) by a Fourier transform with r = |x|. It is the mass parameter m gap that carries the information about the (physical) mass gap of QCD. For example, it leads to the confinementdeconfinement temperature T c (m E gap , m M gap ), where m E,M gap are the chromo-electric (inverse temporal screening length) and chromo-magnetic (inverse spatial screening length) mass gaps at finite temperature. In turn, in particular m E gap is also sensitive to the confinementdeconfinement phase transition, see e.g. [49,50]. In the vanishing temperature limit both masses reduce to m gap , which is thus directly linked to confinement, for more details see [9,10,48].
We emphasize that m 0 may change significantly with only minor or no changes inflicted on m gap for m 0 m gap , and thus the value of m 0 does not influence observables such as T c . This is sketched in Fig. 2 (blue branch) and suggests that m 0 with m 0 > m gap labels different IR solutions of Yang-Mills theory in Landau-type gauges. In turn, the red branch in Fig. 2 with m 0 ≈ m gap → ∞ clearly differentiates between physically different theories labelled by m gap . Evidently, only solutions to functional equations within the former regime may describe Yang-Mills theory, while the latter regime may be interpreted as solutions in massive Yang-Mills theory with physically different masses. Interestingly, the solution seen on the lattice is close to (or at) the boundary between these two regimes; see [42] for more details. Note also that the mass gap in the Yang-Mills branch is the minimal one that can be obtained in the system.
We also note that the scaling for the propagators in (4) can be embedded in an IR hierarchy for all n-point functions, which admits a 1/Q 4 behavior for all diagrams that contribute to the quark-antiquark (qq) four-point function [51,52]. This leads to a linear rise in the qq potential with the distance r in coordinate space already in a single gluon exchange picture. If such a behavior could be shown to be present for r 1/m gap , it would facilitate the access to observables showing direct signatures of confinement. For the decoupling solutions, on the other hand, no order of a diagrammatic expansion of the qq four-point function shows a 1/Q 4 behavior. However, for the whole family of solutions with m 0 m gap , confinement appears at the confinement-deconfinement phase transition with an m 0 -independent T c [9,10]. Moreover, the expectation value of the Polyakov loop, which is the respective order parameter, can be obtained within a resummation of diagrams [53].
The open questions we want to address in this work are therefore: What is the mechanism that generates the IR singularities in Fig. 1? What distinguishes the different decoupling solutions, and is there a preference for one or the other type of solutions?
In the following we show that the emergence of the IR singularity in Fig. 1 is tied to the longitudinal projection of the gluon self-energy, which we call Π(Q 2 ) and which for the decoupling solutions acts as an effective mass term (see Eq. (13) below). In the present work we offer two possible interpretations: In Scenario A, we interpret this term as an artifact of the regularization and/or truncation; in this case we find that only the scaling solution survives, however with an ambiguity in the IR exponent κ. In Scenario B, we consider the term to be dynamical; in this case gauge invariance can only be preserved if there is a massless longitudinal pole either in the ghost-gluon, three-gluon or four-gluon vertex, which drops out from the transverse equations and only serves to eliminate Π(Q 2 ). Here we use the assumption of complete IR dominance of the ghost, and a peculiar finding is the fact that such a pole does appear in the ghost-gluon vertex, but only for the scaling solution.
The paper is organized as follows. In Sec. II we discuss , gluon propagator (middle) and three-gluon vertex (bottom). The gluon self-energy depends on the ghost loop, gluon loop, tadpole, squint and sunset diagrams. The three-gluon vertex contains the ghost triangle, gluon triangle and swordfish diagrams, plus further two-loop terms and diagrams with higher n-point functions (not shown). The ghost-gluon vertex (green) and four-gluon vertex (orange) satisfy their own equations and are external inputs in the system.
the DSEs for the Yang-Mills system in Landau gauge and the different truncations that we employ. In Sec. III we explain Scenario A and its consequences. In Sec. IV we investigate Scenario B and the emergence of massless longitudinal poles in the ghost-gluon vertex. We conclude in Sec. V. To keep the paper self-contained, the explicit diagrams appearing in the DSEs are worked out in detail in Appendix A. Appendices B, C and D provide further details on the renormalization, the numerical procedure and the longitudinal poles in the ghost-gluon vertex. We work in Euclidean conventions throughout the paper.

II. YANG-MILLS DSES
In the present work we aim at a thorough understanding of the mechanisms at work within the mass generation in Yang-Mills theory in Landau gauge QCD in a functional formulation. While quantitatively reliable approximations have been set up in functional approaches, see e.g. [42,64], we shall employ approximations here that carry all the dynamics of the system but are still simple enough to access these dynamics directly.
Accordingly, we solve the coupled DSEs for its lowest n-point functions shown in Fig. 3: the ghost propagator, gluon propagator and the three-gluon vertex. The propagator equations are exact since the gluon DSE also includes the two-loop terms. In the three-gluon vertex DSE we neglect further two-loop terms and diagrams includ-ing vertices without a tree-level counterpart. To keep the discussion transparent, we relegate the explicit formulas to Appendix A and only highlight the most important aspects in the main text.
Throughout this work we keep the ghost-gluon vertex at tree level because this is a good approximation in Landau gauge [13,84]: Here the vertex is finite in the ultraviolet (UV) and does not need to be renormalized, and explicit calculations have shown that the deviations from the tree-level behavior are small [56,64]. The threegluon vertex, on the other hand, is almost completely dominated by its classical tensor structure and has only a mild angular dependence [60]. Thus, in the following we restrict ourselves to the leading tensor and the symmetric limit such that the vertex is represented by a onedimensional function F 3g (Q 2 ). The quantities we compute are thus the ghost dressing G(Q 2 ), the gluon dressing Z(Q 2 ) and the three-gluon vertex dressing F 3g (Q 2 ).
We consider the following truncations of the coupled DSEs in Fig. 3: Setup 1: This is the DSE system solved in [14][15][16], where the three-gluon vertex DSE is bypassed and only the propagator equations are solved, also neglecting the two-loop terms in the gluon DSE. To avoid introducing model input for the three-gluon vertex, one can use the ansatz F 3g (Q 2 ) = G(Q 2 )/Z(Q 2 ), which ensures the correct renormalization and perturbative limit of the vertex, or keep the vertex at tree level with F 3g (Q 2 ) = const. For our purposes, Setup 1 will mainly serve as a reference point since it is well established by now that the three-gluon vertex is suppressed at small momenta and likely has a zero crossing [31, 33, 35, 56-58, 60, 63, 83].

Setup 2:
Here we still neglect the two-loop terms in the gluon DSE but back-couple the three-gluon vertex DSE into the system. For the four-gluon vertex, which now appears as an additional input, we employ its classical tensor multiplied with G(Q 2 ) 2 /Z(Q 2 ), which again ensures the correct renormalization of the vertex and its perturbative limit. We will see below that the error induced by this truncation is at the 10% level.
Setup 3: This corresponds to the full system in Fig. 3 also including the two-loop terms in the gluon DSE. As a consequence, the propagator DSEs are two-loop complete at UV momenta. The remaining inputs are the ghost-gluon and four-gluon vertex where we use the ansätze discussed above. Below we will see that the error in this truncation is at the 3-4% level.
We emphasize that there is no explicit model input in any of these truncations (except for dropping higher vertices and tensor structures) and the only parameter in all cases is the coupling g. Furthermore, the qualitative features found in this work are independent of the truncations and appear in all three setups.
The general form of the gluon DSE is given by where with the transverse projection operator defined in (2). In (7a), the tree-level propagator D µν 0 (Q) is given by Eq. (1) with the replacement Z(Q 2 ) → 1/Z A , where Z A is the gluon renormalization constant. The gluon self-energy Π µν (Q) in (7a) is the sum of the ghost loop, gluon loop, tadpole, squint and sunset diagrams in Fig. 3. The transversality of the gluon self energy in (1) and (7) is a direct consequence of the Slavnov-Taylor identities (STIs) for general covariant gauges, e.g. [85].
In approximations transversality might be lost. Indeed, already in perturbation theory the transversality of (7b) is not present in the single diagrams in Fig. 3, and longitudinal pieces that are related by gauge symmetry cancel between diagrams. This is used in the PT-BFM scheme, together with dimensional regularisation, for a reorganization of classes of diagrams according to transversality; see e.g. [25].
Hence, in order to account for truncation artifacts we allow for a more general form of the gluon self energy, as already done in (1), where we allowed for a momentum-dependent longitudinal dressing L(Q) for the propagator. Note that ∆ 0 (Q 2 ) and Q 2 ∆ L (Q 2 ) decay for large momenta Q 2 Λ 2 , see e.g. [8,42,86]. Inserting Eqs. (1) and (8) in (7a) and comparing coefficients leads to the following equations for the propagator dressing functions G(Q 2 ), Z(Q 2 ) and L(Q 2 ), These relations hold in all linear covariant gauges. The set of DSEs in (9) makes it apparent that the tensor basis used in (8) is over-complete: one of the ∆'s can be reabsorbed in a redefinition of the other two. In the absence of approximations, (7b) holds and we arrive at We introduced the over-complete basis (9) because there are two qualitatively different sources for the ∆'s, which is important for the discussion of the IR behavior: In numerical applications of functional approaches, loop integrals are regularized by a momentum cutoff Λ, where loop momenta k 2 > Λ 2 are dropped. The respective renormalization as well as the related Bogolyubov-Parasiuk-Hepp-Zimmermann (BPHZ) renormalization requires mass counter terms for guaranteeing the cutoff independence and transversality of (7b), see also [87]. Potential remnants are proportional to δ µν and may lead to ∆ 0 . Accordingly, a consistent treatment of these artifacts is of eminent importance for avoiding explicit mass terms in the IR. Any truncation in which tensor structures of vertices are dropped may lead to artifacts, which can be distributed between ∆ T and ∆ L due to the overcompleteness of the basis. While the ∆ T contribution is simply a correction to Π phys , the generation of ∆ L is potentially harmful in the IR.
The generation of ∆ L can be elucidated at the relevant example of the ghost-gluon vertex Γ µ gh . Its complete form has two tensor structures, where p µ is the outgoing ghost momentum and Q µ the incoming gluon momentum. In Landau gauge, the ghost renormalization contant Z Γ = 1 so that the dressing functions A(p 2 , p · Q, Q 2 ) and B(p 2 , p · Q, Q 2 ) measure the deviation from the classical vertex. In our present approximation we set A = 0. Evidently, if we also set B = 0 we would remove a completely longitudinal part of the ghost diagram in the gluon self energy in Fig. 3 adding to ∆ L (Q 2 ). This entails that the cancellation of all longitudinal parts in the sum of diagrams will be absent and the self-energy will have longitudinal parts. Moreover, as B drops for large momenta, so will the longitudinal part inflicted by this approximation. Accordingly, this specific approximation artifact does not affect multiplicative renormalization and in particular does not contribute to the mass counter term. In summary, the regularization with a momentum cutoff as well as truncations may lead to artifacts in the gluon DSE that may complicate the identification of the transverse part, in particular in the IR.
In numerical applications of functional approaches, the over-complete basis in (8) is usually not used; instead one absorbs either ∆ 0 or ∆ L in two other dressing functions. We denote these by with Π(Q 2 ) = Π(Q 2 ) + Π(Q 2 )/Q 2 , which lead to the following self-energy decompositions: Π(Q 2 ) and Π(Q 2 ) correspond to the decomposition which we employ for the explicit formulas in App. A. It facilitates the consistent treatment of the UV renormalization by keeping the δ µν term explicit, and it is convenient because Π(Q 2 ) and Π(Q 2 ) are free of kinematic 1/Q 2 singularities. For the decoupling solutions ∆ 0 (0) becomes constant and ∆ T (Q 2 → 0) has a logarithmic divergence stemming from the ghost loop. Note, however, that within approximations the transverse part Π = ∆ T − ∆ L also contains the purely longitudinal contribution ∆ L , which can come from missing tensor structures in vertices such as the longitudinal part B of the ghost-gluon vertex in (11). A consistent treatment of these terms may be of importance for the generation of the mass gap. The transverse-longitudinal decomposition on the other hand, corresponds to the DSEs for Z(Q 2 ) and L(Q 2 ) in (9). In this case, the transverse part Π(Q 2 ) includes potential artifacts of the UV renormalization, which must be treated consistently to avoid the implicit introduction of an explicit gluon mass.
As an example, we consider a numerical solution of the DSEs projected on the transverse part by contracting (13) or (14) with a transverse projection operator, hence computing the dressing Π(Q 2 ) = ∆ T (Q 2 ) + ∆ 0 (Q 2 )/Q 2 . In Fig. 4 we show exemplary decoupling and scaling solutions obtained in the present work for Π(Q 2 ) and their breakdowns into individual diagrams. The ghost loop contribution is positive, the gluon loop negative and the two-loop terms only have a small effect.
Interestingly, ∆ 0 /Q 2 drops rapidly and already vanishes for values of Q 2 in the mid-momentum region, which entails that the UV renormalization is done consistently. In turn, ∆ 0 /Q 2 contributes to the IR behavior. In fact, in the decoupling case the 1/Q 2 singularity in the inverse gluon dressing and hence the gluon mass gap is solely produced by this term, which is tantamount to an explicit mass. In the scaling case, the term is still present but has a steeper divergence induced by the ghost loop which becomes large. As a consequence, both ∆ 0 /Q 2 and ∆ T diverge with the same power in the IR, which leads to the IR ghost dominance with the scaling behav- Clearly, all these points need to be understood and the consistent definition of the transverse dressing is of eminent importance for the correct description of the generation of the gluon mass gap. Hence, in the following we shall consider two scenarios: Scenario A: We assume ∆ L ≡ 0, i.e., we drop terms that are purely longitudinal. Therefore, we must dispose of ∆ 0 as well to ensure Π = 0. With this assumption we can extract the respective transverse dressing by an appropriate projection. This is discussed in detail in Sec. III.
Scenario B: We assume ∆ L = 0, in which case we need ∆ L = −∆ 0 /Q 2 for consistency, see (10). Then, the transverse dressing is simply obtained by contracting the DSE with T µν Q . This is discussed in detail in Sec. IV and entails a dynamical generation of the mass gap.
We shall see that the generation of the gluon mass gap in Scenario A, in the present approximation, is not dynamical but rather similar to that in massive Yang-Mills theory. In turn, the generation of the gluon mass gap in Scenario B is indeed dynamical. While this entails a clear preference for Scenario B, the discussion in Sec. III will serve to introduce the concepts and provide the technical details, which are the same in both scenarios along with many of the results.
We also note that there are no ambiguities in the practical extraction of any set of two linearly independent dressing functions. The only purely longitudinal part proportional to ∆ L that can appear in our approximations arises from the longitudinal tensor in the ghostgluon vertex, which we treat as an 'extra' term in Sec. IV. Thus, for the explicit formulas in Appendix A without this term, which are based on the decomposition (13), one can identify Π = ∆ T and Π = ∆ 0 .

III. GLUON MASS GAP: SCENARIO A
Scenario A, as defined at the end of the last section, is based on the assumption that purely longitudinal terms are absent, i.e., ∆ L (Q 2 ) = 0 in Eq. (8). Then, Π(Q 2 ) = ∆ 0 (Q 2 ) has to be an artifact either stemming from the hard cutoff employed in the DSEs and/or from the truncation of the equations. Thus, in the full system with a gauge-invariant regulator or with appropriately defined counter terms it would vanish identically. In this way, systematically improving the truncations would improve the precision on Π(Q 2 ) while eventually sending Π(Q 2 ) to zero, so we might as well drop it right away.
A way to interpolate between the two cases with and without this term is to contract the self-energy with the general projection operator with a parameter λ, where the transverse and longitudinal projection operators T µν Q and L µν Q have been defined in (2). Contracting (13) or (14) with (15) leads to where λ = 1 corresponds to the transverse projection and λ = 0 is the Brown-Pennington projection [88].
To begin with, for λ = 0 one cannot find QCD-like decoupling solutions, or even convergent DSE solutions. This is already apparent in the numerical solution in Fig. 4. There, the mass gap and m 0 solely originate from Π. In its absence, the system would have to exhibit a Coulomb-type solution which is not present.
In turn, the scaling solution in principle still exists, but for λ = 0 the analytic IR analysis leads to a 0/0 problem [14,15,43], which is shown in Fig. 5. For a tree-level ghost-gluon vertex, the algebraic condition that determines the IR exponent κ reads For a transverse projection (λ = 1) this yields κ ≈ 0.595. Sending λ → 0 entails κ → 1/2, but for λ = 0 the solution disappears because λ/(1 − 2κ) = 0/0. The above considerations suggest to regularize the limit Π → 0: We replace it by an auxiliary mass term, which is finally sent to zero in a controlled way. Thus, we investigate the situation with a constant mass term Here, g is the coupling and G 0 = G(Q 2 = 0) is the value of the ghost dressing function at the origin. This ensures the correct renormalization as will become clear below. Furthermore, β is a dimensionless mass parameter and µ is the renormalization scale. In this way, the equations explicitly depend on a mass parameter β which is sent to zero in the end. Note that this is formally similar to massive Yang-Mills theory, i.e., the Curci-Ferrari-model which has been studied in a series of works, e.g. [43,71,89]. However, in our case the mass parameter β does not arise from the Lagrangian; it is merely an artifact which must be taken to zero to recover massless Yang-Mills theory. Note that Eq. (18) is also a conceptual simplification since there is no need to deal with quadratic divergences.

A. Renormalization
What can obscure the study of DSEs to some extent is their dependence on the renormalization constants or, equivalently, on the renormalized values of the gluon and ghost dressing functions Z(µ 2 ) = Z µ and G(0) = G 0 . Hence, before we proceed with the explicit numerical solutions in Section III B, we discuss the renormalization of the Yang-Mills system in detail. We note that the following discussion is independent of the two scenarios and equally applies to Scenario B.
Propagators, vertices and the coupling are related to their bare counterparts, denoted by the superscript (B), by the multiplicative renormalization constants, where as a consequence of the STIs. In the Landau gauge the ghost-gluon vertex stays unrenormalized, so we can set Z Γ = 1. Therefore, all renormalization constants can be related to Z A and Z c , We can consistently renormalize the ghost and gluon propagators at different renormalization points: In practice it is convenient to renormalize the gluon dressing function at Q 2 = µ 2 to Z(µ 2 ) = Z µ , while the ghost dressing is fixed at Q 2 = 0 to G(0) = G 0 . This fixes all renormalization constants and the resulting equations depend on g, Z µ and G 0 , see Appendix B for details: Here, Σ is the ghost self-energy, Π = Π + Π/Q 2 the gluon self-energy in Eq. (14), and M are the vertex diagrams. The renormalization constants are dynamically determined by which fixes Z 3g and Z 4g from Eq. (20) so that no further renormalization of the vertices is required.
In practical solutions of the Yang-Mills DSEs one keeps g and Z µ fixed; the value G 0 of the ghost dressing at vanishing momentum then distinguishes the scaling and decoupling solutions. Any finite G 0 corresponds to a decoupling solution and the limit G 0 → ∞ to the scaling solution. This is, however, not completely satisfactory from the viewpoint of renormalization as sketched in Fig. 6. Z µ and G 0 should only multiplicatively renormalize the propagators, which would not lead to physically different solutions; in a double-logarithmic plot, a renormalization only induces vertical shifts in both functions. Similarly, the value of the coupling g should only set the scale, and when plotted over a logarithmic momentum scale it would only induce horizontal shifts in both functions. But how is it then possible to obtain a family of different decoupling solutions?
To this end, let us redefine the renormalized propagators by dividing out their values at the respective renormalization scales, with the same redefinition for all other quantities that renormalize like G and Z, i.e., the renormalization constants and the three-and four gluon vertices This is always possible since the dressing functions are only fixed up to multiplicative renormalization. For the scaling solution, G 0 becomes infinite. However, throughout this work as well as in most of the literature, the scaling case is defined as the limit of the decoupling solutions for G 0 → ∞. Then, all redefinitions are well-defined. By inspection of the DSEs, see Appendix B, one finds that g, Z µ and G 0 no longer appear individually in the equations but only in the combination Henceforth, we call α the coupling parameter, which can take any positive value. Note that α should not be confused with the strong coupling g 2 /(4π); it has the µ running of the gluon propagator and will later be related to the mass parameter. The equations (21)(22) with redefined quantities take the form Here we have extracted the factor α from the redefined self-energies, which do not explicitly depend on it, apart from the two-loop terms. Thus, all renormalization constants have disappeared and only the coupling α remains. With Eq. (26) it becomes clear that it is technically not G 0 that distinguishes the decoupling solutions but α: If g and Z µ are kept fixed, changing G 0 is equivalent to changing α, but equivalently we could have kept any other two variables fixed and changed the third variable. Thus, any finite α produces a decoupling solution and α → ∞ corresponds to the scaling solution.
It is also convenient to remove the explicit dependence on the renormalization scale µ. To do so, we introduce a dimensionless variable x, which is always possible since the momentum scale still carries arbitrary units and only rescales the dressing functions according to Fig. 6. If we perform the same operation for the loop momenta inside the integrals, then the dependence on µ 2 drops out from the equations and moves to the cutoff of the integrals. The mass parameter β disappears from Π(Q 2 ) as well, which in the redefined system becomes In turn, β appears in the subtraction point since Q 2 = µ 2 entails x = 1/β. To give a concrete example, in the redefined system the ghost self-energy from Eq. (A19) in Appendix A reads Here, u = k 2 /(βµ 2 ), u ± = u + x/4 ± √ ux z, and the µ 2 dependence has moved into the cutoff L. Finally, we redefine the ghost dressing function, and every quantity that renormalizes with it, once more by As a consequence, the coupling α no longer appears in front of the self-energies but only in the ghost DSE. Then the renormalized equations are given by with Π(x) = Π(x)+1/x. This is equivalent to the original equations but now the only explicit parameters are α and β. The subtraction point is no longer arbitrary but tied to the mass parameter β. In particular, the massless limit β → 0 corresponds to a subtraction at x → ∞. Here it is also obvious why α → ∞ corresponds to the scaling solution, because α only appears in the ghost DSE and nowhere else, with G(0) = √ α. Note that we could have equally absorbed the coupling in the gluon dressing by setting Z(x) → αZ(x). This would lead to G(0) = 1 and Z(1/β) = α but for 0 < α < ∞ the solutions are always identical up to multiplicative factors. At this point, the scale x in all equations above is still arbitrary and given in internal units. A quantity that is invariant under renormalization is the running coupling which remains finite even for α → ∞. Here the parameter α is absorbed in the ghost dressing by Eq. (31). In dynamical QCD, the scale should be fixed to experiment by setting the value of α(x) at some momentum scale, which sets the scale in GeV units. In Yang-Mills theory, the usual procedure is to set the scale by comparing with quenched lattice QCD. At asymptotically large momenta, the three-and fourgluon vertex dressings satisfy Hence, we can define equivalent 'running couplings' by multiplying α(x) with powers of the renormalization-group invariants At large momenta these couplings are all identical but in the IR they can be different. In our calculations we set the four-gluon vertex to F 4g (x) = G(x) 2 /Z(x), so the second relation in (34) is trivially satisfied. The deviation from the first relation will serve as a measure of the truncation error, which we discuss below.

B. Discussion of solutions
We proceed with the numerical solution of the DSEs (32) in the (α, β) plane, where α is the coupling and β the mass parameter as explained above. The selfenergy contributions for the ghost, gluon and three-gluon vertex DSEs are worked out in Appendix A. For the gluon DSE in Scenario A, we only need to keep the terms contributing to Π(Q 2 ), whereas those for Π(Q 2 ) are replaced by the constant mass term. This also means that the tadpole diagram drops out completely since it contributes to Π(Q 2 ) only. For the same reason we also do not need to worry about quadratic divergences; all self-energy diagrams are only logarithmically divergent and these divergences are removed by the subtraction. Fig. 7 shows the results for the running coupling α(x), the ghost dressing G(x), the gluon dressing Z(x) and the three-gluon vertex dressing F 3g (x) in Setup 3, i.e., with a back-coupled three-gluon vertex and including the twoloop terms in the gluon DSE. The calculations were performed at a fixed value β = 0.001. The variation with β will be discussed in Sec. III C. The different curves correspond to a variation of the coupling parameter α over its full range. The running coupling is invariant under renormalization but the propagator and vertex dressing functions are not, so we plot the bare dressing functions according to Eq. (19) for better visibility, i.e., With increasing but finite α, the ghost dressing function in Fig. 7 rises in the IR and saturates at a finite value. The gluon dressing behaves like Z(x) ∝ x in the IR and develops a maximum at intermediate momenta. As a consequence, also the running coupling x in the IR and develops a maximum, however at a different value of x.
The three-gluon vertex becomes increasingly smaller in the IR, where it eventually crosses zero and diverges logarithmically. The orange bands in Fig. 7 mark the onset of the 'QCD-like' decoupling solutions around α 1, as we discuss below. For α → ∞, the curves eventually converge to the scaling solution, where the dressing functions in the IR behave like In this way, the scaling solution is the envelope of the decoupling solutions and produces a finite running coupling α(x) which becomes constant in the IR. Thus, the emergence of a family of decoupling solutions is simply a consequence of varying the coupling parameter α, and the scaling solution is the limiting value for α → ∞.
The coupled DSEs provide us with an internal way to quantify the truncation error. In practice the equations do not converge unless we modify the renormalization constant Z 3g in Eq. (20) by a parameter c, This is equivalent to an additional renormalization of the three-gluon vertex in Eq. (21) of the form Note that this additional renormalization would not be necessary if the STIs were preserved. Since we are working in different truncations defined by the Setups 1, 2 and 3, the STIs are no longer satisfied, but the effect can be compensated by introducing the parameter c = 1. Note however that the STIs constitute a further system of functional relations, and in non-trivial nonperturbative approximations such as the present truncations of the vertex expansions it is impossible to satisfy all sets of functional relations simultaneously; see e.g. [7,8,39,42] for respective discussions. From Eqs. (33)(34)(35), the parameter c 2 translates to the ratio α 3g (x)/α(x) evaluated in the UV. This reflects the consistency of the running couplings in the UV, which is known to be important for the quantitative accuracy of the solutions. It turns out that for each truncation there is a value c max where the anomalous dimensions of the ghost and gluon propagators reach their physical values. In each iteration step we employ fits of the form (40) where the coefficients and exponents are free fit parameters. Thus, they are not forced to reproduce the anomalous dimensions of Yang-Mills theory given by The rightmost plot in Fig. 8 shows the maximum value of the gluon dressing function Z max = max Z (B) (x) as a function of c. Close to c max , there is a strong curvature which even appears to tend to infinity for the scaling solution. The horizontal orange band corresponds to the region of lattice results for Z max , see below. The extrapolated values are compatible with this band, although the extrapolation induces a large uncertainty.
The results in Setups 1 and 2 are qualitatively similar to those in  Without any truncation one must have c = 1, so the value of c max allows us to quantify the truncation error. With a tree-level three-gluon vertex it is about 40%, which reduces to ∼ 10% if the three-gluon vertex is backcoupled and 3 − 4% if also the two-loop terms are included. With further systematic improvements such as a back-coupling of the ghost-gluon and four-gluon vertices, the error should thus become even smaller. Indeed, then the respective c, or rather the ratio of all different couplings, approaches unity and the deviation is at the sub-percent level, see [42,64].
Close to the respective value c max , the dressing functions in Fig. 7 are very similar in all setups. In doublelogarithmic plots like those for α(x), G(x) and F 3g (x) (recall that asinh F grows logarithmically for large |F |) they are hardly distinguishable, except that in Setup 1 the three-gluon vertex stays at tree-level. The quantity that is most sensitive to the truncations is the bump of the gluon dressing function whose height defines Z max . In Setup 2 the bump is slightly larger and narrower than in Setup 3 because the two-loop terms in the gluon selfenergy are positive (cf. Fig. 4) and therefore flatten Z(x).
In any case, these observations imply that the gluon dressing function may well have only reached a fraction of its true height. The largest possible value c = 0.95 we can reach in Setup 3 corresponds to the results in Fig. 7. The same observation holds true, although less pronounced, for the ghost dressing at intermediate momenta. This makes it difficult to compare with lattice results since we should match the results at c = c max but the functions can change substantially in its vicinity. In addition, the lattice data are only determined up to scale setting (which amounts to horizontal shifts of the curves on a logarithmic scale) and multiplicative renormalization. Finally, we must identify the value of α that should be compared to the lattice-type decoupling solution.
In practice we construct simple parametrizations for the lattice data in Ref. [34], The fits in (42) implement the one-loop running with γ gh , γ gl from Eq. (41), together with G lat (0) = 1 and Z lat (x → 0) ∝ x. They are determined up to normalization and rescaling and must be matched to the UV running of the DSE results. The behavior of the UV coefficient a gh in Eq. (40) is analogous to Fig. 8. The curves converge approximately to the same value at c max independently of the coupling α, from where we can match the ghost dressing with the lattice. Because this changes G lat (0) and our G(x) satisfies G(0) = √ α, with this we identify α ∼ 2.5 as the value corresponding to the lattice decoupling solutions. Matching the UV running of the gluon together with the condition Z(1/β) = 1 then fixes the gluon dressing function and the scale. From the resulting curves in Fig. 9 one can see that there is a substantial gap between the lattice data and the DSE results obtained with c = 0.95. However, we find a similar behavior when we apply naive cubic extrapolations like those in Fig. 8 to G(x) and Z(x) over the whole range in x: The extrapolated ghost dressing at c = 0.97 reproduces the lattice data, and the extrapolated gluon dressing matches the height of the lattice data but its bump is somewhat shifted.
To check the self-consistency of our comparison, we also solved the ghost DSE for G(x) in Eq. (32), with Σ(x) from Eq. (30), using a fixed gluon input. In this case the standalone ghost DSE converges without problems, since the gluon propagator is just an input function. For that purpose, we generated a family of curves for Z(x) which interpolate between the DSE and lattice results (shown in the right panel of Fig. 9). The corresponding results for G(x) from the ghost DSE are shown in the left panel. One can see that the DSE solution with the lattice gluon input reproduces the lattice ghost without difficulties. The ghost DSE only depends on α. Solving it for different values of α, we find that the quantity x |G(x) − G lat (x)| is minimized for α ≈ 2.5 in agreement with our findings above. This also confirms that the only missing ingredient of the ghost DSE, the dressing of the ghost-gluon vertex, can only have a minor effect.
The problem of not being able to reach the physical value c max in the present truncation appears to be related to the back-coupling of the three-gluon vertex and, by extension, to the four-gluon vertex which appears in its respective DSE. Our results are stable when changing the number of grid points or using different grids, but achieving convergence close to c max becomes increasingly difficult even with relaxation and Newton methods (cf. App. C). We easily obtain large and narrow bumps for Z(x) when we employ model ansätze for F 3g (x), but this comes at the expense of additional parameters and scales. This would obscure the discussion of the results in view of the gluon mass gap, and we refrained from doing so. We can also reduce the gap between the lattice and DSE results for Z(x) by identifying the lattice decoupling solution with DSE solutions at different values of α, but in those cases we can no longer simultaneously describe the ghost dressing function.
Concerning the four-gluon vertex ansatz, we employed F 4g (x) = G(x) 2 /Z(x) but when we use instead F 4g (x) = Z 4g = const. the results are very similar. However, in the vicinity of c max the shape of F 4g (x) especially in the mid-momentum region appears to play an increasingly prominent role in (de-)stabilizing the solutions. In Fig. 3 one can see that the four-gluon vertex appears in the gluon DSE only in the sunset diagram, which turns out to be negligibly small, so it mainly enters indirectly through the swordfish diagram in the three-gluon vertex. The general four-gluon vertex involves many tensors, including three momentum-independent ones which may lead to additional effects [64,69]. For functional computations aiming at quantitative precision in Yang-Mills theory, which also implement a dynamically back-coupled four-gluon vertex, we refer to Refs. [42,64]; e.g. in [64] the deviation from the STI was quantified to be on the sub-percent level.
In summary, our analysis implies that the lattice results can be matched to the DSE results for a particular value of the coupling parameter α. In the next section we will see that this generalizes to a particular line of constant physics in the (α, β) plane.

C. Lines of constant physics
The results discussed so far were obtained at a fixed value of the mass parameter β = 0.001. The remaining question in Scenario A is how the respective results and findings change with β. For modest changes with larger β, we find that the plots in Fig. 7 are essentially unchanged except they may be shifted horizontally in log x. The overall scale is arbitrary and we are free to rescale the solutions, and after rescaling they are almost identical. Thus, in this domain the mass parameter β rescales the solutions but does not change any physics. This is in contradistinction to a variation of the coupling parameter α, which distinguishes physically different solutions.
For quantifying these observations, we show the lines of constant physics in the (α, β) plane in Fig. 10. For their determination we extracted the maximum of the running coupling α(x) for each parameter set (α, β). Then, we rescaled x such that the maximum always appears at the same point x = x 0 . Accordingly, a given line is defined by the same value α(x 0 ). This amounts to interpolating the dressing functions at fixed value of β = 0.001 and finding the value of α for a given α(x 0 ). Along each line, all functions are identical up to rescalings, which can be seen in Fig. 11: every curve therein is a superposition of 30 curves along a given line and these curves are practically indistinguishable.
In Fig. 10 one can see that for sufficiently large β the lines of constant physics become indeed horizontal. On the other hand, if they were horizontal for all values of β this would be incompatible with physical expectations: In the massless theory with β = 0, the coupling α should not produce physically different solutions but only rescale the system. In other words, in the massless limit the line of constant physics should be the vertical axis in Fig. 10. These different forms can only be reconciled with each other if the lines bend towards the origin in the (α, β) plane. Indeed, this is what we find for very small values of β, as is visible in Fig. 10. The behavior in Figs. 10 and 11 is completely analogous in Setups 1 and 2 and thus appears to be a general feature. In practice we cannot reach β = 0 for the reasons discussed in Appendix C. 1/β is the subtraction point in the gluon DSE, and the shape of the gluon self-energy defines a calculable window for β. For increasing values of α or c, the gluon dressing becomes increasingly narrow and this window shrinks. The smallest value we can reach is typically β ∼ 10 −6 . The results in the left panel of Fig. 10 were obtained with c = 0.92; in the right panel one can see that the effect becomes more pronounced when c is increased towards its practical limit c = 0.95, i.e., the curves become steeper. However, for c 0.92 we can no longer cover the full range in β.
The interpretation of our results is sketched in Fig. 12. The two parameters (α, β) in the Yang-Mills system form a combination λ 0 (α, β) that acts along the lines of constant physics and only rescales the system without changing any physics. The second combination λ 1 (α, β) is orthogonal to these lines and distinguishes the physically different solutions, i.e., only λ 1 is an active parameter. In other words, in the presence of a mass term the cou- pling parameter α no longer rescales the system but distinguishes the physically different decoupling solutions, with the scaling solution as their endpoint.
In the limit β → 0, the shape of the lines implies that λ 0 (α, 0) = α, i.e., the coupling α only rescales the solutions while not changing the physics. However, this means that the scaling solution obtained at β = 0 and α → ∞ is identical to the scaling solution at β = 0 with arbitrary α, and it is the only solution that remains in the massless limit. This implies that in the present Scenario A the only physical solution is the scaling solution, as the mass parameter β is merely an artifact that must be sent to zero. The scaling solution is therefore a parameterfree, intrinsic property of the system. Moreover, since the scaling solution as well as the whole family of decoupling solutions implement confinement via the gluon mass gap [9,10], this implies confinement in the Landau gauge.
While this seems theoretically appealing, there is a conceptual problem that can be traced back to the determination of the IR exponent κ discussed around Eq. (17). We have replaced Π(Q 2 ) by a constant mass term of the form (18), which is 'removed' in the end by sending β → 0. As β only appears in the term β/x, we could absorb it into the scale and thus sending β → 0 is equivalent to renormalizing the equations at x = 1/β → ∞. This entails, however, that the mass term 1/x cannot truly be removed as shown in Fig. 13, which is the decomposition of the gluon self-energy Π(x) in Scenario A. We emphasize that its analogue in Fig. 4 corresponds to Scenario B, which we discuss below. For fixed β and α → ∞ the ghost loop becomes large, but even for the scaling solution at α → ∞ the term 1/x is present and dominates the IR behavior. Had we not absorbed β into the scale, for β → 0 the same would happen but all curves would be pushed to the left on a logarithmic scale, so that the rise with 1/x in the IR eventually drops out of the numerical grid and is shifted to log x → −∞. Therefore, the effect of the mass term is never truly switched off.
The same rescaling effect explains the apparent discrepancy between the solutions for Z(x) at α → 0 from Eqs. (27) and (32), where the first set of equations entails Z(x) → 1 and the other Z(x) → 1/(1 + 1/x − β). In the former case, the structure effects are also pushed to log x → −∞, so the renormalized theory never becomes a free massless theory for α → 0. Because the 1/x term drives the IR properties, in this case the scaling solution at x → 0 with has an IR exponent κ = 1/2. This is also what we found analytically in Fig. (5) when sending the parameter λ to zero. The ambiguity for λ = 0 reappears in the numerical solution in the following way: Since the mass term in Eq. (18) is arbitary as long as we send β → 0 in the end, we are free to replace it with any other ansatz of the form where κ is a free parameter, as long as we remove this term again in the end by sending β → 0. After the redefinitions from Sec. III A this becomes which replaces the 1/x term in the DSEs (32). Its effect is again negligible in the mid-and large-momentum regions, so it does not change the behavior of the solutions, but as long as κ ≥ 1/2 it still dominates the IR and leads to a scaling solution with IR exponent κ. We checked this numerically and found that by dialling κ we can indeed generate any scaling solution we want, where the lines of constant physics behave qualitatively in the same way as before.
In other words, in Scenario A the limit β → 0 is not unique, which is a manifestation of the 0/0 problem in the analytic IR analysis and leads to a 'family of scaling solutions' depending on the arbitrary value of the IR exponent κ. Because this ambiguity does not seem very appealing, in the next section we revisit the initial problem and discuss a possible alternative scenario.

IV. GLUON MASS GAP: SCENARIO B
The second interpretation, which we call Scenario B, is that ∆ 0 (Q 2 ) is not an artifact but indeed a dynamical feature of the equations. To study it, we return to Eq. (12) with both ∆ 0 (Q 2 ) and ∆ L (Q 2 ) in the system. Then the transverse projection to arrive at Π(x) is achieved with λ = 1 in Eq. (16). This is the transverse system that is commonly solved in DSE and fRG studies. However, in view of the previous discussion, it begs a few questions which we investigate below: the appearance of quadratic divergences in Sec. IV A and, after a brief discussion of the solutions, the issue of gauge consistency in Sec. IV C, i.e., the vanishing ofΠ.

A. Quadratic divergences
The first and more practical issue is that a dynamical ∆ 0 term potentially comes with quadratic divergences. While ∆ 0 vanishes in perturbation theory using dimensional regularization, a hard cutoff interferes with gauge invariance so that ∆ 0 has an intrinsic 'artifact' admixture that must be removed. Once the quadratic divergences are subtracted in a fully consistent way, the remainder (if it is non-zero) must be a dynamical, nonperturbative effect.
Different methods have been employed to remove quadratic divergences in the sum Π = ∆ T + ∆ 0 /Q 2 , see [90] for an overview. One is an explicit subtraction at the level of the integrands; in that case one analyzes their UV behavior and subtracts the terms producing the quadratic divergences. This has been employed in what we call Setup 1 [15] but it becomes cumbersome when the three-gluon vertex is back-coupled dynamically or two-loop terms are included. Another method is to subtract the quadratic divergences numerically by fitting Π to 1/Q 2 terms [91]. A consistent combination of both is typically done in the fRG, where the consistency follows from the (trivial) RG consistency of the fRG setup [8,42].
Having individual access to ∆ T and ∆ 0 , a simpler method is to subtract ∆ 0 directly, i.e.
Since ∆ 0 (Q 2 0 ) is a constant, the arbitrariness in the choice of Q 2 0 is absorbed by adding another arbitrary constant proportional to βµ 2 , again with a prefactor that ensures the correct renormalization. This leads to a similar form of the equations as in Scenario A, where we added a constant mass term with mass parameter β. After the redefinitions leading to Eq. (32), the self-energy becomes with σ(x) given by This procedure is equivalent to the one in Ref. [64], where a term ∝ 1/x is subtracted from the self-energy. In this way, the only difference between the Scenarios A and B is the term σ(x). Note that we can no longer change the power of the 1/x term like we did in Eq. (45) in Scenario A, because the subtraction must be done in its numerator. As a consequence, the IR exponent κ is no longer arbitrary but a dynamical result with a definite value depending on the truncation. Because σ(x)/x is only active in the IR, x 0 should not be too large to ensure numerical stability; in practice we choose x 0 = 1. The arbitrariness of x 0 is then absorbed in the change of β which determines the subtraction point x = 1/β.

B. Discussion of solutions
The numerical results in Scenario B do not require a separate discussion, because the only change in the equations is the addition of the term σ(x) in Eq. (47a), which only affects the IR region. Thus, in the mid-and highmomentum regions, and also at low momenta as long as α is not too large, all previous statements remain intact. The shape of the solutions in Fig. 7, the behavior of the UV exponents in Fig. 8, and the values of c max for the different setups all remain approximately the same in Scenario B. The only noticeable difference appears in the IR, where σ(x) is active. Fig. 14 shows the gluon self-energy contributions in this case. These are the same plots as in Fig. 4 except now we separate the 1/x and σ(x)/x contributions explicitly. For the decoupling solutions, σ(0) is a constant. The remaining self-energy terms in ∆ T (x) are either constant or diverge at best logarithmically in the IR (see App. A), so that Π(x) and thus Z(x) −1 are dominated by the terms ∝ 1/x for x → 0. These terms generate a mass gap such that the (dimensionless) gluon propagator at the origin reads After reinstating dimensions through x = Q 2 /(βµ 2 ), the r.h.s. is divided by the (arbitrary) factor βµ 2 , from where one can extract the gluon mass scale according to Eq. (3). When approaching the scaling solution for α → ∞, σ(x → 0) grows and eventually diverges with a power x 1−2κ stemming from its ghost loop contribution. At the same time the ghost loop contribution to ∆ T (x) also diverges with x −2κ in the scaling limit. Taken together, both ∆ T (x) and σ(x)/x contribute with the same power to the IR divergence of Π(x) ∝ x −2κ , which for κ > 1/2 beats the divergence from the mass term 1/x as is visible in the right panel of Fig. 14.
As a result, the dressing functions in the IR scale with The IR exponent κ can be read off from any of the dressing functions in the scaling limit. This is shown in the left panel of Fig. 15 for the ghost dressing, where the line 1/x κ is the envelope of the decoupling solutions. In practice we determine κ from the ghost self-energy by fitting it to Σ(x) = a + bx κ for x → 0, where a, b and κ are free fit parameters depending on {α, β, c}, i.e., also for the decoupling solutions. From Eq. (32) we deduce that the ghost dressing function in the IR is well approximated by so that scaling at x → 0 is achieved for α → ∞ whereas for each decoupling solution G(0) = const. For sufficiently large α we may read off the scaling exponent κ ≈ 0.58, which is compatible with contemporary DSE [63] and fRG results [42]. Note however that the approximate value undershoots the analytic scaling κ = 0.595 that is obtained in the present approximation for α → ∞. In the right panel of Fig. 15 one can see that κ is almost independent of c, i.e., the fact that we cannot exhaust the full range up to c max ≈ 0.96 . . . 0.97 is irrelevant for its determination. When we solve the DSEs in the (α, β) plane we find the same behavior as in Scenario A. One can identify lines of constant physics, which for β → 0 bend towards the origin as shown in Fig. 16. The interpretation is thus the same: One combination of α and β rescales the solutions and another one distinguishes the decoupling solutions, with the scaling solution as their endpoint. The crucial difference to Scenario A, however, is the fact that β is no longer an artificial mass parameter that must be removed in the end by sending β → 0, because now it arises from the subtraction of quadratic divergences. Therefore, any value for β is equally acceptable, which also seems more appealing considering the conceptual problems with the limit β = 0 encountered in Scenario A. As a consequence, there is no longer a criterion that discriminates between the different solutions. Instead of an ambiguity in the IR exponent κ, one is left with the ambiguity between the scaling and decoupling solutions, which for β = 0 are distinguished by the parameter α. This raises the question: Can one find another criterion that allows for a selection of solutions?

C. Longitudinal singularities
The fundamental question in Scenario B is howΠ(Q 2 ) can vanish, given that the dynamical term ∆ 0 (Q 2 ) = 0 now contributes to it. As we already mentioned, this implies a non-vanishing ∆ L (Q 2 ) for the consistency relation ∆ L = −∆ 0 /Q 2 to hold. To this end, let us re-examine the assumptions we made. The crucial (implicit) assumption in our approximation so far is the use of the ghostgluon vertex at tree-level. When introducing Scenario B in Section II, we already emphasized that the ghost-gluon vertex has two tensor structures, cf. Eq. (11): The momentum p µ is the outgoing ghost momentum and Q µ is the incoming gluon momentum. The corresponding dressing functions are A(p 2 , p·Q, Q 2 ) and B(p 2 , p·Q, Q 2 ). Accordingly, the tensor structure missing in the present approximation produces to a solely longitudinal contribution to the inverse gluon propagator which belongs to the ∆ L part. The ghost-gluon vertex is UV-finite in Landau gauge and does not need to be renormalized. The function A is known to be small [56,64], which is why a treelevel vertex is a good approximation. However, little is known about the function B apart from the limit p µ = Q µ where the incoming ghost momentum vanishes, in which case the vertex becomes bare [84] and thus one has B(p 2 , p 2 , p 2 ) = −A(p 2 , p 2 , p 2 ).
The Lorentz tensor Q µ attached to B is purely longitudinal and vanishes whenever the gluon leg is contracted with a transverse gluon propagator in Landau gauge, i.e., for all internal gluon lines. As such, it does not contribute to the ghost self-energy in Fig. 3, but it does contribute to the ghost loop in the gluon DSE. Because the term is longitudinal, it drops out from the transverse projection Π = ∆ T + ∆ 0 /Q 2 in Eq. (12). Hence, the full transverse dressing is not affected by the non-classical dressing B in the ghost-gluon vertex (modulo small effects from A = 0). However, the longitudinal projection picks up an extra term and generalizes to where ∆ 0 is the sum of all Π-like self-energy contributions listed in Appendix A (up to A = 0 effects). The longitudinal dressing ∆ L can be worked out in analogy to Eq. (A21): where k ± = k ± Q/2. The resulting DSEs read However, Eq. (52) allows us to eliminate the completely longitudinal term Π(Q 2 ) as required by the STIs. For Π = 0 we must have For decoupling solutions, ∆ 0 (0) is a constant. Thus, Eq. (55) can only hold for all Q 2 if the longitudinal ghostgluon dressing B diverges like 1/Q 2 for Q 2 → 0: the ghost-gluon vertex must exhibit a longitudinal massless pole. In this case Π vanishes, and we arrive at the purely transverse self energy with Π = ∆ T + ∆ 0 /Q 2 . The equation for Z(Q 2 ) is unchanged and B only serves to eliminate Π(Q 2 ).

Complete infrared dominance of the ghost
Note that ∆ L (Q 2 ) in Eq. (53) only constitutes the part of the longitudinal dressing that originates in the ghost loop. However, while other diagrams in principle may also contribute to ∆ L (Q 2 ), the dominant contribution in the IR limit Q 2 → 0 is given by (53). In particular, for the scaling solution all the other terms rise with lower powers of 1/Q 2 . In turn, this complete IR dominance of the ghost part for the scaling solution is only successively weakened for the family of decoupling solutions.
For completeness we also consider the sub-leading contributions from the gluon loops. In principle, a longitudinal tensor in the ghost-gluon vertex also propagates to the three-gluon vertex through the ghost triangle. However, we find that in the symmetric limit for the vertex such a term only gives a contribution to one out of three possible tensors. This contribution drops out from the gluon DSE because in the Landau gauge it is contracted with internal transverse gluon lines, see the discussion around Eqs. (A56) and (A70) in the appendix. Therefore, in the truncations we consider herein such a longitudinal term can only affect the ghost loop in the gluon self-energy, i.e., it can only arise from the ghost-gluon vertex. This argument also further corroborates the IR dominance of the ghost contribution in the gluon self energy for decoupling solutions that are sufficiently close to the scaling limit. Based on this argument we proceed with the assumption of complete IR dominance of the ghost contributions.
If there are indeed longitudinal poles in the ghost-gluon vertex, they should dynamically emerge in its DSE which is shown in Fig. 17. While we cannot make any statement about the ghost-gluon four-point function in the last diagram therein (for a corresponding discussion, see Ref. [92]), it is clear that the longitudinal tensor must drop out for the internal vertices since they are contracted with transverse gluon lines. Thus, the B term can only survive in the upper vertex in the second diagram (marked in yellow). This yields an inhomogeneous Bethe-Salpeter equation (BSE) for B with the structural form (see Appendix D for details) Here, KB represents the second (Abelian-like) diagram in the r.h.s. of Fig. 17 and the inhomogeneous term B 0 is the sum of the remaining loop diagrams; the tree-level term does not contribute to B. Eq. (57) will produce a singularity in B at some value i if the kernel satisfies K = 1, i.e. if any of its eigenvalues satisfy λ i (Q 2 = −m 2 i ) = 1. In that case, the numerator defines a Bethe-Salpeter amplitude ϕ via If one is only interested in the pole locations, one can equivalently solve the homogeneous BSE Kϕ i = λ i ϕ i , which does not require knowledge of B 0 and depends on the Abelian-like diagram in Fig. 17 only. The typical eigenvalue spectrum of a BSE kernel is sketched in Fig. 18. The eigenvalues λ i (Q 2 ) decrease with Q 2 , which is also true in our case because the ghost propagators in the loop depend on the external gluon momentum and fall off at large Q 2 . For Q 2 > 0, the eigenvalues must be smooth functions of Q 2 because the integrand does not have spacelike singularities. Each intersection with λ i (Q 2 = −m 2 i ) = 1 corresponds to a singularity of B at Q 2 = −m 2 i , which may appear at timelike values Q 2 < 0. The largest eigenvalue λ 0 (the 'ground state') corresponds to the pole closest to the origin in Q 2 . In particular, if λ 0 satisfies λ 0 (Q 2 = 0) = 1, then B must have a massless pole.
We can then solve the homogeneous BSE directly at Q µ = 0 and calculate λ 0 (0). Setting A = 0 in Eq. (51) to be consistent with our truncation of the DSEs, the resulting equation becomes very simple, Here we employed the same redefinitions and rescaling that led to Eq. (32) for the coupled DSEs, i.e., x is the dimensionless variable corresponding to p 2 . The kernel only depends on the ghost and gluon dressing functions G and Z.
The left panel of Fig. 19 shows the largest eigenvalue λ 0 (Q 2 = 0) of K as a function of the coupling parameter α. For α → 0, it is clear from Eq. (59) that the eigenvalue must be zero because the product ZG 2 contains an intrinsic factor α. On the other hand, the limiting value is λ 0 (0) = 1 because λ 0 (0) > 1 would produce a pole for spacelike values Q 2 > 0 and thus a tachyonic singularity. One can see in Fig. 19 that for all decoupling solutions with finite α one has λ 0 = 1. However, as α rises and approaches the scaling solution at α → ∞, the eigenvalue approaches λ 0 → 1. Thus, we find that the scaling solution indeed has a longitudinal massless singularity in the ghost-gluon vertex but the decoupling solutions do not. As a consequence, for the decoupling case the condition (55), which ensures Π(Q 2 ) = 0 and therefore gauge invariance, cannot hold. In conclusion, in Scenario B with the assumption of complete infrared dominance of the ghost only the scaling solution survives.
The results in Fig. 19 have been obtained within the Setup 3, but we find the same behavior in Setups 1 and 2. The result is also independent of the parameter c, which is consistent with the observation that the scaling properties are independent of c (cf. Fig. 15). In the analogous plot for Scenario A, λ 0 does not reach the value 1 but this is also not necessary because the condition (55) does not arise.
We also note that a necessary requirement for a massless singularity in the ghost-gluon vertex is the fact that the vertex in the present formulation is not completely ghost-antighost symmetric; see e.g. [47,90] for more general discussions. If that were the case, KB could not produce a massless pole strong enough to ensure Eq. (55) (see Appendix D for details). Another remark is that the ghost-gluon vertex DSE in Fig. 17 is the so-called 'c DSE', where the vertex on the left in each diagram is bare. An equivalent form is the 'A DSE', where all top vertices are bare (see e.g. Fig. 2 in [63]). In the latter case one would not arrive at a BSE for the B term and the longitudinal poles would need to come from higher n-point functions such as the four-ghost or two-ghost-two-gluon vertices in the t channel [92].

General case
The results of the last section beg the question whether the present setup or an (implicit) IR completion of the Landau gauge only admits the scaling solution. Indeed, the complete IR dominance of the ghost assumed above is only valid in the scaling limit for α → ∞. Accordingly, one should interpret the above result as a corroboration of a gauge-consistent scaling solution and not as an exclusion of gauge-consistent decoupling solutions. For the latter solutions there is no IR suppression of the gluonic contributions with powers of Q 2 , and one may have to also include the three-gluon vertex in the analysis. The coupled BSE system of ghost-gluon and threegluon vertex has been studied intensively for lattice-type decoupling solutions in the PT-BFM scheme [21,22,[24][25][26][27][28]. As discussed earlier in Sec. III B, in our present setup the lattice solution roughly corresponds to α Lat = 2.5 at β = 0.001, which is shown by the vertical dashed line in the left of Fig. 19. In particular, in [27], within the approximations used there, an almost complete dominance of the gluonic contributions has been found. This may entail that for α Lat the BSE for the ghost-gluon vertex is only fully consistent if also including the three-gluon term in Fig. 17. Then, the BSE in (57) generalizes to where K gh-gl is the contribution from the first diagram in Fig. 17 that we considered so far. In turn, K 3gl stands for the second diagram in Fig. 17 with the three-gluon vertex. The ghost-gluon BSE in (60) has to be augmented with that for the three-gluon vertex and has to be solved simultaneously for general α. While being of eminent importance, a detailed analysis is deferred to future work. Here we simply put forward a possible scenario that encompasses the present consistent scaling solution with complete IR ghost dominance and the lattice-type decoupling solution in [27] with relative three-gluon dominance. In combination, the present findings and that in [27] suggest that the massless singularity is present in the system for α α Lat . Moreover, within a crude linear analysis the value of λ 0 (Q 2 = 0) for the lattice-type decoupling solution with α ≈ α Lat , cf. left panel in Fig. 19, gives a rather small estimate for the contribution of K gh-gl to the required full eigenvalue λ 0 (Q 2 = 0) = 1, which suggests the dominance of the gluonic contribution. Such a scenario is fully compatible with the findings in [27].
In turn, in massive Yang-Mills theory with a large explicit mass m 0 ≈ m gap → ∞ we do not expect such a massless singularity to be present. Moreover, as discussed earlier, α Lat corresponding to the lattice solution is slightly above the onset of the decoupling solutions around α 0 ∼ 1 highlighted in Figs. 7 and 8 and thus the massive Yang-Mills regime is defined by α α 0 . The same statement within a quantitative approximation has been shown in [42].
In combination, this makes us speculate about the exciting possibility to distinguish the possible QCD-type solutions with α α 0 from that in massive Yang-Mills theory with α α 0 by the existence or absence of the massless pole in the longitudinal sector. This presence of this singularity is indicated by and the scenario described above is sketched in the right panel in Fig. 19. Whether it can be validated within more sophisticated truncations remains to be seen.

Other scenarios
It has also been speculated that Π(Q 2 ) could be nonzero because the longitudinal DSE in (54) depends on the combination ξ Π(Q 2 )/Q 2 . Since ξ = 0 in Landau gauge, one would obtain L(Q 2 ) = 1 even if Π(Q 2 ) = 0. However, even if Π were non-zero in Landau gauge, L = 1 must still hold in all other linear covariant gauges and thus Π would still need to vanish for ξ = 0, in which case it is a non-analytic function of the gauge parameter ξ.
This scenario is in direct contradiction with the STI: While it is seemingly natural in terms of the gluon propagator, the STI requires the gluon self-energy (13) to be transverse, and hence Π = 0.

V. SUMMARY AND CONCLUSIONS
In this work we studied the coupled Dyson-Schwinger equations (DSEs) for the ghost propagator, gluon propagator and three-gluon vertex in the Yang-Mills sector of QCD in Landau gauge using different truncations. We addressed several questions related to the origin of mass generation.
We clarified the role of renormalization and showed that the parameter that distinguishes the decoupling so-lutions is the coupling α, which no longer rescales the system in the presence of a mass term. Such a dynamical mass term arises from the longitudinal projection of the gluon self-energy, which we called Π(Q 2 ) and which in principle appears to break gauge invariance. To this end, we offered two scenarios to remedy the problem.
In Scenario A, we considered the mass term to be an artifact of the hard cutoff and/or truncation, so we replaced it by a constant mass parameter β and studied the limit β → 0. We investigated the lines of constant physics in the (α, β) plane and found that one combination of α and β rescales the solutions while the other distinguishes the physically different decoupling solutions, with the scaling solution as their endpoint for β → 0. While the limit β → 0 cannot be reached in practice, we nonetheless conclude that in this case only the scaling solution survives, however with an ambiguity in the infrared scaling exponent κ.
In Scenario B, we did not modify the DSEs by hand and the dynamical mass term remains in the system. As a consequence, one needs to subtract quadratic divergences which introduces again an arbitrary mass parameter β. However, in this case the longitudinal projection of the gluon self-energy can only vanish if either of the vertices that enter in the gluon DSE has longitudinal massless poles. In the three truncations we studied, this leaves the ghost-gluon vertex as the only candidate. Using the assumption of complete infrared dominance of the ghost, we found that the vertex indeed has such a pole, which is only present for the scaling solution.
Indeed, this solution is the only one which exhibits complete infrared dominance of the ghost. Interestingly, the massless pole is not present within this complete ghost dominance setup for any decoupling solutions. However, for lattice-like decoupling solutions massless poles have been found within the PT-BFM scheme. In contradistinction to the scaling case, they are dominated by the gluonic contributions. The combined observations allowed us to put forward a scenario that encompasses the present consistent scaling solution with complete infrared ghost dominance and the lattice-type decoupling solution with three-gluon vertex dominance, for more details see Section IV C 2. A detailed numerical study of this scenario will be presented elsewhere. In this appendix we collect the expressions for the ghost propagator, gluon propagator and three-gluon vertex DSEs. The propagator DSEs in Fig. 20 read where the full propagators depend on the ghost and gluon dressing functions G(Q 2 ) and Z(Q 2 ): Their tree-level expressions (with subscript 0) follow from the replacements G(Q 2 ) → 1/Z c and Z(Q 2 ) → 1/Z A , where Z c and Z A are the ghost and gluon renormalization constants, respectively. The transverse and longitudinal projectors are given by T µν Q = δ µν − Q µ Q ν /Q 2 and L µν Q = Q µ Q ν /Q 2 and ξ is the gauge parameter. Writing the ghost self-energy as Σ G (Q) = −Q 2 Σ(Q 2 ) and decomposing the gluon self-energy as in Eq. (13), the DSEs take the form where the functions Π(Q 2 ) and Π(Q 2 ) are obtained from Eqs. (15)(16) using λ = 0. According to the decomposition (13), we will refer to Π(Q 2 ) as the transverse part and Π(Q 2 ) as the gauge part in what follows. Note however that for Π(Q 2 ) = 0 the transverse part Π differs from the transverse projection Π = Π + Π/Q 2 . Every self-energy diagram contains a tree-level vertex, whose Feynman rules are collected in Fig. 21 and Table I. The full three-gluon vertex depends on 14 tensors [93]; here we restrict ourselves to the classical tensor in Eq. (A12), which amounts to the replacement In addition, we assume that the dressing function only depends on the symmetric variable: Concerning the four-gluon vertex, it is convenient to rewrite the tree-level expression (A13) as where the Lorentz and color tensors are given by Because (τ 3 ) ab,cd = f ade f cbe is linearly dependent due to the Jacobi identity (τ 3 = τ 1 − τ 2 ), Eqs. (A9) and (A13) are identical.
The full four-gluon vertex depends on 136 Lorentz tensors [69] and five color structures. Here we keep again only the tree-level tensor with a general dressing function F 4g obtained by the replacement Also in this case we assume that the dressing function only depends on the symmetric variable:  In the following we work out the self-energy diagrams in Fig. 20 explicitly. We work in Landau gauge (ξ = 0) where the ghost-gluon vertex is UV-finite and thus we set Z Γ = 1.
In hyperspherical variables, the Lorentz-invariant integral measure k = d 4 k/(2π) 4 takes the form where Λ is the hard cutoff. Since the only Lorentz invariants appearing in the single-loop diagrams are Q 2 , k 2 and ω = k · Q = √ k 2 Q 2 z, the integrations over the variables y and ψ become trivial. This corresponds to the frame where The internal loop momenta in Fig. 20 are k µ ± = k µ ±Q µ /2 such that k 2 ± = k 2 + Q 2 /4 ± ω.
The ghost self-energy becomes where the sign in front comes from the DSE (the selfenergy appears with a minus sign). With Eqs. (A3-A4) and f acd f bcd = N c δ ab , the resulting expression is Note that the self-energy is negative, as it should be according to the DSE (A5): the ghost dressing function is enhanced compared to the tree-level expression and the inverse dressing function is suppressed. For Q 2 = 0, the self-energy integral becomes constant both in the scaling and decoupling case.
The ghost loop in the gluon DSE is given by In principle there are two signs in front, one from the minus of the self-energy and the other from the closed fermion loop, which cancel out. The sign of −f acd f bcd = −N c δ ab finally cancels the sign coming from the vertices. By applying Eqs. (15)(16) with λ = 0 we extract the transverse and gauge parts: In principle there are also terms ∼ z in the bracket, but these vanish after integration since the integral over z is symmetric. Note that the term with (1 − 4z 2 ) does not produce a quadratic divergence at large k 2 or singularities at small Q 2 : in that case one has k 2 ± ≈ k 2 and the integration over z yields Thus, only the gauge part Π gh has a quadratic divergence whereas the transverse part Π gh is only logarithmically divergent. Its contribution to the gluon DSE is positive, so it has the tendency to reduce the gluon dressing function compared to its tree-level expression. For Q 2 → 0 and the scaling solution, both Π gh and Π gh /Q 2 diverge with the same power (Q 2 ) −2κ , whereas for the decoupling solutions Π gh diverges logarithmically but Π gh remains finite.
The gluon loop in the gluon DSE reads where the minus in front comes from the sign of the selfenergy and the factor 1/2 from the symmetrization. Taking again the projections and abbreviating Z ± = Z(k 2 ± ) and F 3g = F 3g (k 2 − , k 2 + , Q 2 ), we arrive at with the kernels and Once again, quadratic divergences only appear in the gauge part Π gl whereas the transverse part Π gl is only logarithmically divergent due to the factor (1 − 4z 2 ). It is negative and has the tendency to increase the gluon dressing function. For Q 2 → 0, Π gl and Π gl become constant for both scaling and decoupling solutions.
The tadpole diagram is given by which is very simple to work out since only the tree-level four-gluon vertex appears: The transverse part vanishes due to the factor (1 − 4z 2 ), cf. Eq. (A22), which is not only true in Landau gauge but also for a general gauge parameter ξ. Thus, the tadpole only gives a contribution to the gauge part, which is a constant and quadratically divergent. Because in Scenario A the gauge part is replaced by a constant, and in Scenario B a constant is subtracted from the gauge part by Eq. (47a), the tadpole drops out from all equations.

Gluon DSE: Two-loop terms
The two-loop terms in the gluon DSE are small compared to the one-loop terms (see Fig. 4). Nevertheless, they complete the DSE and ensure that the ghost and gluon propagators are two-loop exact in the UV. In addition to the loop momentum k from Eq. (A17), we write the second loop momentum l as The propagator momenta are then given by l ± = l ± k/2 and k = k − Q as shown in Fig. 22. In this case the innermost integration is trivial, so that there are in total five integrations instead of two. The integral measure is This would lead to a substantially higher computational demand, but with the simplifications for the three-and four-gluon vertices explained around Eqs. (A8) and (A15) the integrands of the two-loop terms simplify substantially. To this end, we define the variables where a hat denotes a normalized momentum, with the inverse relations Then we can rewrite the two innermost integrations as Below we will see that the integrands only depend on Ω but not on w, so that the w integration becomes trivial. In addition, with some rearrangements of the integrals it turns out that the squint diagram requires the same computational demand as the one-loop terms and the sunset diagram involves only one additional integration. The squint diagram gives the dominant contribution to the two-loop terms. It reads where the prefactor −1/2 comes from the minus sign in the DSE and the symmetrization, and the transverse projectors from the gluon propagators were defined below Eq. (A4). We absorb all Lorentz-invariant dressing functions into a quantity I sq (l 2 , k 2 , Q 2 , z, Ω) defined by where the three-gluon vertices depend on the symmetric variables The squint diagram then becomes where the factor 27/2 is the color trace and the remaining Lorentz kernel K µν sq is just a kinematic expression. With the decomposition (15)(16) we can split it into its Lorentzinvariant components: × K sq (Q 2 , k 2 , l 2 , z, z , y) K sq (Q 2 , k 2 , l 2 , z, z , y) .
Using Eq. (A33), and because I sq does not depend on the variable w, we can integrate out w to arrive at The kernels K sq and K sq are given by with the coefficients Even though one must integrate over k 2 , l 2 , z and Ω, the product of the second, third and fourth line in Eq. (A39) can be written as and The functions A and A only require looping over three variables (Q 2 , k 2 and z) and the same is true for B(k 2 ) which does not depend on Q 2 , hence the loops go over k 2 , l 2 and Ω. Therefore, calculating the squint diagram is computationally no more expensive than the one-loop diagrams, at least within our approximations for the threegluon vertex. Once again, Π sq (Q 2 ) is only logarithmically divergent and Π sq (Q 2 ) diverges quadratically.
The sunset diagram is almost negligible compared to the squint diagram, which dominates the contribution from the two-loop terms by far. It is given by where the prefactor −1/6 comes from the sign in the DSE and the symmetrization factor. We employ the same strategy as before and absorb the dressing functions in a quantity I sun (l 2 , k 2 , Q 2 , z, Ω) defined by where the four-gluon vertex dressing depends on the symmetric variable The sunset diagram can then be written as with the color trace 9/2 and a kernel K µν sun that is purely kinematic. The resulting dressing functions are Π sun (Q 2 ) Π sun (Q 2 ) = − 1 6 g 4 Z 4g 9 2 k l I sun (l 2 , k 2 , Q 2 , z, Ω) × K sun (Q 2 , k 2 , l 2 , z, z , y) K sun (Q 2 , k 2 , l 2 , z, z , y) .
Once again, I sun does not depend on the variable w, which can be integrated over to arrive at In this case the kernels K sun and K sun are more complicated, but they are still even functions of Ω which involve at most quartic powers (1, Ω 2 and Ω 4 ). Thus we may expand them as follows: and likewise for K sun . The product of the second, third and fourth line in Eq. (A48) then takes the form where A n (Q 2 , k 2 , l 2 ) The coefficients c n read explicitly: Also in this case, the potentially dangerous terms proportional to k 2 /Q 2 do not produce quadratic divergences since they come with factors (1 − 4z 2 ), and the terms proportional to k/Q are multiplied by z so that their integrals vanish for k 2 → ∞ or l 2 → ∞. The coefficients entering in Π sun (Q 2 ), on the other hand, do lead to quadratic divergences: We note that the {A n , A n } do not depend on Ω and the B n do not depend on z and Q 2 . Thus, we still saved two integrations such that the sunset diagram only requires one additional integration compared to the one-loop and squint diagrams.

Three-gluon vertex DSE
Because we back-couple the three-gluon vertex into the propagator DSEs, we must also work out its own DSE. It is given by and the corresponding diagrams are shown in Fig. 23.
Here we neglected further two-loop diagrams as well as those containing higher n-point functions (such as the two-ghost-two-gluon vertex). We also restrict ourselves to the symmetric limit, where p 2 1 = p 2 2 = p 2 3 = Q 2 , and to the classical tensor structure from Eq. (A13). As shown in Ref. [60], this is a good approximation since the angular dependence is mild and the effect from higher tensors are small. A symmetrization of the diagrams is not necessary in the symmetric limit where all six permutations are identical, hence we contract two of the three swordfish diagrams from Fig. 3 into one (the factor 2 is implicit).
The full three-gluon vertex in the symmetric limit has the general form and depends on three fully antisymmetric tensors, where q i = p j − p k and {i, j, k} is an even permutation of {1, 2, 3}. K 1 is the classical tensor, whereas the remaining ones carry higher momentum powers and thus their dressing functions are subleading. Moreover, K 3 vanishes upon transverse projection and drops out from the coupled DSEs in Landau gauge where it is always contracted with two transverse gluons in the loops. Eq. (A55) can be verified as follows. In general, the three-gluon vertex has 14 possible tensors, which are arranged in Table V of Ref. [60] in terms of singlets, antisinglets and doublets under the permutation group S 3 . Since the full vertex is Bose-symmetric and the color structure f abc antisymmetric, the Lorentz part must be antisymmetric as well. The three tensors K 1 , K 2 and K 3 are already antisymmetric (they correspond to A (ψ 1 ), A (ψ 2 ) and A (ψ 4 ) in the table), whereas the remaining ones would have to be combined with the momentum doublet D = {a, s} to construct antisymmetric tensors, where the variables a and s are given by However, these variables vanish in the symmetric limit and thus D = 0, so that in the symmetric limit only the three tensors above survive. In the following we work out M µνρ diag (p 1 , p 2 , p 3 ), where 'diag' stands for the ghost triangle, gluon triangle and the swordfish diagrams in Fig. 23. In practice we expand each diagram in the full basis We also note that if we had implemented the full ghostgluon vertex from Eq. (51), then by the projection (A59) the longitudinal B term (which may have massless poles) contributes only to M (3) gh (Q 2 ), which is the dressing function of the tensor K 3 in Eq. (A56). However, K 3 is longitudinal and drops out from the DSEs in Landau gauge because it is contracted with two internal transverse gluon lines in each loop diagram where the threegluon vertex appears. Therefore, the longitudinal poles from the ghost-gluon vertex do not couple into the threegluon vertex, at least not in the symmetric limit. For the same reason K 3 would drop out from the gluon DSE and the ghost-gluon vertex DSE in Fig. 17. In our setup, longitudinal poles can thus only emerge from the ghostgluon vertex.
For the gluon triangle the color factor is the same apart from a minus sign, cf. Eq. (A12), and there is no prefactor in front of the integral. This yields where the indices under the bar are contracted with the transverse projectors from the gluon propagators with respect to the first two momentum arguments of each vertex. In this case the tree-level component becomes where the kernel K gl is given by Here we abbreviated γ = 1 − z 2 and the coefficients are The gluon triangle contribution to the three-gluon vertex is usually positive and small.
The two swordfish diagrams in Fig. 23 turn out to be the leading loop contributions to the three-gluon vertex DSE, apart from the ghost triangle which dominates the IR but is otherwise small (see Fig. 8 in [60] and Fig. 24 in Ref. [63] for similar results; in our Setups 2 and 3 we find the gluon triangle to be even more suppressed). Each of them has a symmetry factor 1 2 and the first diagram picks up a factor 2 because it counts twice in the symmetrization. Using f acd f bcd = N c δ ab and then the combination of color factors from the three-and four-gluon vertex results in the four-gluon vertex combination Γ 1 + Γ 2 /2 from Eq. (A9), which is equivalent to an effective vertex Γ µνρσ 4g (p 1 , p 2 , p 3 , p 4 ) = F 4g (δ µρ δ νσ − δ νρ δ µσ ) (A76) to be multiplied with 3 2 g 2 N c . The first swordfish diagram contains the tree-level four-gluon vertex (F 4g → Z 4 ) whereas in the second diagram the full vertex F 4g (x 4 ) depends on the symmetric variable x 4 . Together with all prefactors, we arrive at where the barred indices denote again transverse projection. Here the tree-level component becomes where the kernel is given by (1 − z 2 ) 4k 2 (3 + y 2 ) + Q 2 (5 − y 2 ) . (A79) Observe that the y dependence inside the integrand is only carried by K sf , so we can integrate it out and use instead K sf = 1 2 dy K sf = k 2 18 (20k 2 + 7Q 2 )(1 − z 2 ) . (A80) K sf is positive and thus the contribution from the swordfish diagrams is usually negative.
Solving the DSEs in Fig. 3 is numerically quite nontrivial because the dressing functions can diverge both in the IR and UV, change signs, and they generally emerge from competing contributions which can vary substantially over many orders of magnitude.
The prime example is the gluon DSE in Eq. (32), where one must ensure that its r.h.s. is positive for all x in each iteration step, since a zero crossing would lead to a spacelike pole in the gluon dressing function Z(x) and thus a tachyonic pole in the gluon propagator, i.e., The self-energy Π(x) is sketched in Fig. 24. It diverges in the IR; with increasing x it eventually becomes negative and develops a minimum before at large x it approaches zero from below. This leads to the condition which defines an interval for the subtraction point x sub = 1/β. Thus, β can neither become too small nor too large.
In addition, the interval can differ in the iteration steps since Π(x) can take different intermediate values before convergence is reached. Ideally the starting guess for Π(x) should not be too different from the converged solution to prevent the function from doing wild jumps during the iteration. In practice we always scan the full α range from α → 0 up to α → ∞, where we extrapolate the converged solutions at the previous values of α as the starting guess for the next α. For small α convergence is easy to achieve, whereas with increasing α the bump of the gluon dressing which defines the shape in Fig. 24 becomes increasingly narrow and this limits the window in β. For intermediate values of α the typical intervals are 10 −6 β 0.5. Because convergence can be very slow, we employ a Newton method whose implementation in the Yang-Mills system was first described in Ref. [94]. To this end, we first map the momentum interval x ∈ (0, ∞) to an interval ω ∈ (−1, 1) in such a way that we can adjust the IR and UV cutoffs ω min and ω max and adjust the density of grid points in the IR, UV and mid-momentum regions separately. Next, we expand the dressing functions G(x), Z(x) and F 3g (x) in polynomials, ln G(x) −1 = n g n P n (ω) , ln Z(x) −1 = n z n P n (ω) , asinh F 3g (x) = n f n P n (ω) , where the logarithm and asinh ensure that the functions do not vary overly strongly over the momentum range and thus a polynomial expansion converges rapidly (n max = 64 is usually sufficient). In practice we employ Legendre polynomials for P n (ω) although the choice does not matter. The inverse relations are g n = dω Ω(ω) P n (ω) ln G(x) −1 , z n = dω Ω(ω) P n (ω) ln Z(x) −1 , where Ω(ω) are the polynomial weights according to the orthogonality relation dω Ω(ω) P m (ω) P n (ω) = δ mn .
The DSEs in Eq. (32) then take the form x = r(x), where the vector x = {g n , z n , f n , Z c , Z A } contains the moments of the dressing functions together with the renormalization constants, and r(x) denotes the righthand side of the DSEs. The 'natural iteration method' means that if x s is the result of the current iteration, one calculates x new = r(x s ), implements the vector x new as the new start guess in the r.h.s. and repeats, hoping that the system eventually converges. The convergence criterion is to minimize the function The Newton method tries to find a better guess for the minimum from the linear approximation This finds the direction δx of the largest variation, following a tangent along the minimization function until it crosses zero. We optimize the procedure using backtracking, where we test different solutions x new (λ) = x s + λ δx; we then employ the value of λ that minimizes F (x) 2 along that direction as the starting guess for the next iteration and repeat until the minimum is found. In practice we still solve the DSEs for the dressing functions G(x), Z(x) and F 3g (x) and only insert the Newton step between iterations: After one iteration step in Eq. (32), we convert the solutions to the vector F (x s ) and the Jacobian matrix J, determine x new (λ) and convert it back to G(x), Z(x), F 3g (x), Z c and Z A to be used in the next iteration. Each λ amounts to another iteration, however with the same J. Therefore, the main additional complication is the computation of J.
As an example, consider the ghost DSE in (32) whose self-energy is given in Eq. (A19). After redefinitions and rescaling, it becomes where denotes the (dimensionless) four-momentum integration, N = 4πN c , and u ± = u + x/4 ± √ ux z. The ghost part of the vector F (x) in Eq. (C6) is where the bar denotes the outcome on the r.h.s. of the DSE. The corresponding entries in the Jacobian are together with The ingredients are the derivatives of the ghost selfenergy, where the derivative with respect to f j vanishes because Σ(x) does not depend on the three-gluon vertex. FromḠ(x) −1 = Z c + Σ(x) we have ∂/∂Z cḠ (x) −1 = 1, and the derivative with respect to Z A vanishes because Z A does not explicitly appear in the ghost DSE. From Eq. (C3) we have and therefore ∂ ∂g j Σ(x) = N K Σ P j (ω + ) , where the mapping x ↔ ω implies the same mapping for u ± ↔ ω ± . Thus, the derivatives of the self-energies are structurally similar to the self-energies except for the appearance of the polynomials. The remaining entries of the matrix J are constructed along the same lines. Employing these techniques, the solution of the coupled DSEs for a given parameter set (α, β) typically takes 5 . . . 10 minutes on a single CPU, with 10 . . . 20 iterations until convergence is reached. and b 1 = b(k 2 , k · Q, Q 2 ) are the dressing functions for the top vertex, and a 2 is attached to the vertex on the right. Because the internal gluons are transverse, the contribution from b 2 drops out. Writing the momenta as where the inhomogeneity [. . . ] is finite at Q 2 = 0. As discussed below Eq. (58), one can solve a homogeneous BSE for the Bethe-Salpeter amplitude ϕ, which is the residue at the pole; if that equation has a solution (i.e., if some eigenvalue becomes 1) at a certain value of Q 2 , then there must be a corresponding pole in b.
We can solve the homogeneous BSE directly at Q µ = 0, in which case the function I becomes Since the ghost and gluon dressing functions therein were obtained from the Yang-Mills DSEs with a tree-level ghost-gluon vertex, we set a 2 = 1 for consistency. The resulting homogeneous BSE at Q 2 = 0 is ϕ(q 2 ) = g 2 N c 2 q 2 k G(k 2 ) 2 k 2 Z(l 2 ) l 4 (1−Ω 2 ) ϕ(k 2 ) . (D13) Note that this implies ϕ(q 2 → 0) ∝ q 2 , which is consistent with b ∝ q 2 /Q 2 being a dimensionless function. For large q 2 we find ϕ(q 2 ) ∝ 1/q 2 and the integral is convergent. After employing the redefinitions leading to Eqs. (32), the BSE becomes identical to Eq. (59) in the main text.
Returning to the question of a ghost-antighost symmetric ghost-gluon vertex, in that case the most general form of the vertex is Eq. (D1) with b = (q · Q) b , where the tree-level vertex corresponds to a = 1 and b = 0. The tree-level contribution to the equation for b therefore vanishes and the integrand of the BSE (D13) picks up a factor k · Q/q · Q = kz /(qz), which strongly suppresses the kernel. In particular, we find that the maximum value of the eigenvalue λ 0 for the scaling solution (α → ∞ in Fig. 19) becomes λ 0 ≈ 0.13, so there is no pole at Q 2 = 0. Even if there was a pole, the divergence of the function b ∝ q · Q/Q 2 would not be strong enough to ensure the validity of Eq. (55). The emergence of longitudinal poles in the ghost-gluon vertex can thus be related to the lack of the ghost-antighost symmetry.