{'fraction_non_alphanumeric': 0.04242835562755873, 'fraction_numerical': 0.01269785668936226, 'mean_word_length': 4.438504583879033, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A covariant variational approach to Yang-Mills Theory

We investigate the low-order Green's functions of SU(N) Yang-Mills theory in Landau gauge, using a covariant variational principle based on the effective action formalism. Employing an approximation to the Faddeev-Popov determinant established previously in the Hamiltonian approach in Coulomb gauge leads to a closed set of integral equations for the ghost and gluon propagator. We carry out the renormalization and the infrared analysis of this system of equations. Finally, we solve the renormalized system numerically and compare with lattice results and other functional approaches.


I. INTRODUCTION
The low-order Green's functions of Yang-Mills theory have been the focus of many investigations, both in the continuum and on the lattice. Functional methods such as the functional renormalization group (FRG) flow equations [1] and Dyson-Schwinger equations (DSE) [2] initially concentrated their studies on the case of covariant gauges, since this is the natural choice in a covariant setup, and allows to use BRST symmetry and the ensuing Slavnov-Taylor identities to guide and improve the analysis. Moreover, the Kuga-Ojima criterion [3,4] argues that there should be a direct link between the deep infrared behavior of the gluon and ghost propagator, and the issue of colour confinement. Most investigations initially found an infrared vanishing, scaling type of solution for the gluon propagator, which is, however, at odds with high-precision lattice simulations [5][6][7][8]. It was only later realized that infrared finite so-called decoupling solutions could also be obtained in the functional approach under certain circumstances [1]. In fact, an infrared finite gluon propagator had been found before in Refs. [9][10][11][12][13][14][15][16][17][18][19][20] By contrast, the so-called variational approach to the Hamiltonian formulation of Yang-Mills theory investigated in [21][22][23][24][25] has been formulated in Coulomb gauge. This non-covariant condition is advantagous in the Hamiltonian approach since it allows for an explicit resolution of Gauß's law. In the Hamiltonian approach of [21][22][23][24][25] the Schrödinger equation for the Yang-Mills vaccum wave functional is approximately solved by the variational principle, minimizing the energy with suitable trial Ansätze for the vacuum wave functional. Compared to functional methods based on the Lagrangian formulation of quantum field theory, like the FRG [1] and DSE [2], the variational principle of the Hamiltonian approach has the benefit that the relative size of the energy density controls the quality of the approximation made, i.e. of the trial wave functional used. By enlarging the space of trial states the description can be improved. The Hamiltonian approach in Coulomb gauge provides direct access to the so-called Coulomb potential between static charges and gives relatively simple explanations for confinement [26][27][28][29][30] and other low-energy phenomena [31]. Furthermore it yields propagators which are in good agreement with lattice calculations [32][33][34][35] and the Gribov-Zwanziger confinement scenario [36,37].
Unfortunately it cannot be generalized directly to covariant gauges and make contact with the alternative studies mentioned above.
In the present paper, we will demonstrate that a conceptually similar variational principle can be established in covariant gauges if one relies on the effective action instead of the energy. We will present a variational approach to quantum field theory which is based on the minimization of the effective action and apply it to Yang-Mills theory in Landau gauge. This approach will result in a system of integral equations for the low-order Green's functions, which we solve using the techniques borrowed from the Coulomb gauge calculations mentioned above. We will also argue that our method gives optimal results for the propagators of the theory, while the extension to realistic vertices presumably requires to go beyond the Gaussian Ansatz using Dyson-Schwinger equation techniques [38]; see also Ref. [39] for an alternative approach.
The paper is organized as follows: In the next section, we present the general variational principle for the effective action and explain in section III how it can be applied to Yang-Mills theory in covariant gauges. After presenting our ansatz for the trial path integral measure, we discuss an approximation to the Faddeev-Popov determinant which was introduced in Ref. [24] in the context of the Hamiltonian approach in Coulomb gauge and which facilitates the treatment of the ghost sector. In section IV we apply the variation principle to the effective action of Yang-Mills theory in Landau gauge and derive a closed set of integral equations for the ghost and gluon propagator. The rather cavalier method to renormalize these equations is put on a solid ground in section V by relating it to the traditional introduction of counter terms. After carrying out the infrared analysis of our renormalized system of equations, we present in section VI their numerical solution and compare with recent high-precision lattice data. Finally, we give a short summary and end with an outlook on further developments.

II. THE VARIATIONAL PRINCIPLE
Below we recall the variational principle for the effective action in quantum field theory. This material mainly serves to fix some of our notation and introduces various concepts used lateron.
We start with a theory for a quantum field φ defined by an action S(φ) in d = n + 1 dimensional Euclidean space time. Expectation values of operators are computed with the normalized measure dµ 0 (φ) = Z −1 exp − −1 S(φ) dφ (1) Here dφ is the flat (translationally invariant) measure in field space, and an implicit regularization is understood. The obvious analogy with statistical mechanics can now be exploited to define a variation method: The Gibbs-like measure eq. (1) is the unique probability measure in field space which minimizes the free energy (or rather, free action), Here . . . µ ≡ dµ(φ) . . . is the expectation value in the trial measure µ, and the normalization is such that 1 µ = dµ(φ) = 1. If the measure is written in Radon-Nikodym form dµ(φ) = dφ ρ(φ) with a suitable density ρ, the quantity is the entropy which measures the accessible field space for quantum fluctuations. The value of the minimal free action from eq. (2) attained for the Gibbs measure (1) is F (µ 0 ) = − ln Z. Notice that Planck's constant really plays a role similar to the temperature in statistical mechanics, i.e. it controls the balance between the classical action S and the fluctuation entropy W (µ).
Traditionally, one would now take the minimising measure eq. (1), or some approximate minimum in a restricted measure space, and compute the Schwinger functions or rather the generating functional W (j) of their connected part. From this, the effective action Γ follows by Legendre transformation.
Alternatively, the effective action can also be characterized directly by a variation principle: To this end, we go back to eq. (2) and perform the minimization in two steps: First, we define a constrained free action by restricting some operator Ω(φ) to a classical value ω, and second, we minimise this action under variation of all trial probability measures, which yields the effective action for the operator Ω, It depends, of course, on the choice of the operator Ω, and is a functional of the prescribed value ω both explicitly (via the constraint) and implicitly (via the ω-dependence of the solving measure). In most cases, the constraint is chosen as the vev of the quantum field itself, so that The variation principle eq. (2) is now equivalent to Γ(ϕ) ! = min, which is a problem in classical field theory that can be considered as solved. We now have two distinct definitions of the effective action: (i) Functional: Define the Gibbs measure as solution of eq. (2), compute the Green functions, or their generating functional, from eq. (4), and finally construct the effective action as generating functional of 1PI correlators by means of a Legendre transformation.
(ii) Linear response: Compute the effective action directly from the variation principle eq. (6) and obtain the 1PI proper n-point functions as derivatives at ϕ = 0.
It is not hard to see that the two descriptions agree, cf. appendix A, but this identity only holds for the exact solution of the variation problem, when no restrictions are placed on the trial measures µ. This is rarely ever the case. In practice, a viable variation scheme will have to restrict the space of trial measures {µ} to those candidates for which the expectation values in eq. (2) or eq. (6) can be computed explicitly. In such a restricted variation, the effective action determined from (i) and (ii) will differ. In appendix B, we sketch the two approaches for a simple φ 4 theory with Gaussian trial measures, and also compare with the variational solution of the Schrödinger equation in the Hamilton formalism.
As with all variational methods, it is not a priori clear which of the two formulations above will give the better approximation to the true system, although general arguments [40] indicate that the linear response approach (ii) is of higher order in the difference between true and approximate minimal measure µ 0 . (However, higher order does not always mean higher accuracy.) In any case, the representation eq. (6) where A and F A are the differential forms for the gauge connection and its field strength, respectively, and g is the bare coupling strength. As usual, Feynman gauge (ξ = 1) simplifies the Lorentz structure of the propagator, while Landau gauge (ξ = 0) yields transversal gluons that can be compared directly to lattice investigations. The prefactors in the action and the inner product with J (0) = 1. It is helpful to interpret J (A) as the weight of the gauge orbit through A, i.e. the canonical volume form on this orbit in field space.

A. Modified variational principle
The appearance of the Faddeev-Popov determinant in the measure eq. (8) requires some modification of the basic variation principle eq. (2), because the latter only holds for measures of the Gibbs form eq. (1). One obvious solution is to transfer the Faddeev-Popov determinant into the action, Then the variation principle based onS takes the standard form F (µ) Here, ρ is the deviation of the trial measure from the flat measure (excluding the Faddeev-Popov determinant), dµ = dA ρ(A). The rhs. of eq. (11) suggests to rewrite the variation principle by redefining the entropy, Here,ρ is now the deviation from the standard measure including the Faddeev-Popov determinant, which is the natural metric on the space of gauge orbits, The redefined entropy eq. (12) coincides with the usual notion of the relative entropy of the trial volume form dA ρ(A) compared to the standard weight dAJ (A) on the space of gauge orbits.
As a consequence, the general variational approach sketched above remains valid for YM theory in covariant gauges, if we only replace the entropy W by the relative entropy W .

B. Gaussian trial measure
In the next step, we have to choose a class of suitable trial measures which is simple enough to allow for the necessary integrals to be performed, but still captures the essential physics of the system. For this purpose, we note Refs. [23][24][25]. The first condition implies that the trial action should be (close to) Gaussian, while the enhancement at the horizon is controlled by the volume form on the gauge orbit: Since the natural volume form is J (A), replacing it in the trial measure by J 1−2α with α ≥ 0 will enhance the weight of near-horizon configurations by a relative factor J (A) −2α ≫ 1.
We will thus attempt a variational approach based on trial measures of the form ( = 1) where N α is the overall normalisation. 1 Note that for α = 0 the Gaussian represents the relative weightρ(A) while for α = 1 2 it gives the full weight ρ(A), cf. eq. (12). We will treat α as a variational parameter, although we shall find below that the exact value of α is immaterial, at least up to two loop order in a formal loop counting scheme introduced in the next subsection.
The measure (14) is unconstrained and thus appropriate for the functional approach discussed 1 The normalisation factor Nα will, in general, depend on the variational kernel ω. In the case α = 1 2 , for instance, we have N 1 2 = det ω/(2π) The normalisation in eq. (14) is such that the bare gluon propagator reduces for g → 0 (and hence J → 1) to i.e. ω(k) → k 2 . This is also the UV limit of ω(k) due to asymptotic freedom. Here, the Lorentz structure is given by the inverse for Landau gauge (ξ = 0) this becomes the transversal projector.
A possible caveat against the trial measure eq. (14) is that it does not respect the BRST symmetry of the full theory, nor any of the identities that follow from it. This is, in a sense, unavoidable in a variational ansatz, because the simplest non-topological action with full BRST symmetry is already the full YM theory, so that any truncation will necessarily break the Slavnov-Taylor identities to a certain extent. For Landau gauge, however, recent lattice calculations favour a soft BRST breaking massive gluon propagator in the deep infrared, and such decoupling solutions were also found under certain assumptions within functional approaches [9][10][11][12][13][14][15][16][17][18]. Since a dynamical mass generation is one of the main virtues of variational methods, the ansatz eq. (14) seems therefore justified for covariant gauges.

C. Curvature approximation
The Faddeev-Popov determinant eq. (10) and its expectation value in the trial measure (14) cannot be computed in closed form. In the following, we will adopt an approximation that has been shown to be correct up to two-loop order in the energy functional within the variational Hamiltonian approach [24] in Coulomb gauge: Since (ln J ) and δ ln J /δA both vanish at A = 0, with a symmetric kernel K[A] that may depend arbitrarily on the gauge connection A. As a consequence, where each digit stands for the combination of colour, Lorentz and spacetime indices, and repeated indices are summed or integrated over. If we now formally introduce a loop counting parameter in the exponent of the trial measure (14), we find from eqs. (19) and (20) ln where the expectation value is taken with the trial measure, eq. (14), and the dots indicate higher orders in the loop counting scheme. To leading loop order, we can therefore replace the kernel by the average curvature on the gauge orbit 2 , In total, this approximation then yields the formula [24] ln which differs, on average, from the exact expression (19) only by a higher loop effect. 3 However, the omitted higher loop terms need not be negative definite and we cannot guarantee that our approximate effective action is always a strict upper bound to the true effective action. We will refer to eq. (22) in the following as the curvature approximation.
The salient point of this approximation is now that the curvature eq. (21) is easier to compute than the expectation value of the full Faddeev-Popov determinant. As shown in section IV C, the curvature can be related to the ghost propagator (although we do not explicitly introduce ghosts), which in turn can be evaluated from a Dyson equation that involves the kernel ω and the full ghost-gluon vertex. In the rainbow-ladder approximation, this full vertex is replaced by the bare one and χ becomes a well-defined function of the ghost propagator and the variation parameters ω and α; for further details cf. section IV C. This approximation is further supported by the fact that the ghost-gluon vertex in Landau gauge is not renormalized [41] and shows little dressing in lattice simulations [42], i.e. the radiative corrections to the vertex tend to cancel.
We can now use the curvature approximation eq. (22) directly in the trial measure (14). From global colour and Lorentz symmetry of our variational ansatz (in the functional formulation without an external classical field A), the expectation value in the definition eq. (21) entails that the curvature has the same simple Lorentz structure as the variation kernel where χ(k) is known as the scalar curvature. The trial measure eq. (14) thus depends on the curvature and the parameter α only in the combination This has the same effect as putting α = 1 2 and replacing ω →ω in eq. (14), i.e. the Faddeev-Popov determinant drops out from our variational ansatz within the curvature approximation which then becomes a simple Gaussian with kernelω(k), and hence results in the gluon propagator (cf. eq. (16)), This observation greatly simplifies the computation of expectation values, since we then have Wick's theorem at our disposal. It must be stressed, however, that the variation is still with respect to ω(k), notω(k), because the curvature χ(k) is, in principle, a dependent quantity.
So far, we have mainly discussed the unconstrained measure for the functional approach, when no external classical field A is prescribed. As discussed in appendix B, this is entirely sufficient to investigate the propagators, since these agree in the functional and linear response approach.
Although we will not study the non-trivial vertex corrections arising in the linear response formulation in any detail, we still want to show at least how the variation problem can be set up in this case: First, we have to adjust the trial measures introduced above to comply with the constraint A = A imposed by the classical field. This can always be achieved by shifting the gauge field A µ → A µ − A µ in the (full) density ρ(A) of the trial measure. 4 In the present case, the density ρ(A) is Gaussian after applying the curvature approximation to eq. (14), so that the final form of our trial probability measure for the linear response approach is . (27) It should be emphasized again that the optimal kernelω A determined from eq. (27) will depend implicitly on the classical field A, which is externally prescribed and thus arbitrary. As a consequence, we can no longer assume the simple colour and Lorentz structure (15) valid in the functional approach; instead we would have to deal with the full matrix gap equation in position space, cf. appendix C.

IV. THE EFFECTIVE ACTION OF YANG-MILLS THEORY
We are now in a position to determine the effective action of Yang-Mills theory from the variational principle using the trial probability measure (27).

A. The free action
To evaluate the free action of the trial measure (27), we first expand the classical Yang-Mills Lagrangian including the gauge fixing term, The relevant correlators in the measure eq. (27) can easily be worked out using Wick's theorem, Next, we have to insert this into S gf and combine it with the relative entropy (cf. eq. (32) below) to obtain the free action F (ω, A). The resulting expressions are, however, rather complicated because the kernelω in the linear response approach does not have the colour and Lorentz symmetry eq. (15). On the other hand, this complication is unnecessary: as argued in appendix B, it is sufficient to use the functional approach with A = 0 and the symmetric kernel eq. (15), as long as we are only interested in the propagators of the theory. For completeness, we present the full expression for the free action F (ω, A) as well as the ensuing gap equation in appendix C, but we do not investigate the resulting vertex corrections in more detail.
For the remainder of this paper, we therefore concentrate on the propagators, i.e. we set A = 0 and use translational invariance to transform to momentum space as in eq. (15). The global colour symmetryω ab ∼ δ ab combined with the anti-symmetry of the structure constants f abc then allows to perform all colour traces by means of the the SU(N ) relations The result for the average classical action is Here, V d is the spacetime volume and the numerical factors are where the symbol d ξ in the first line means d for all ξ = 0, and (d − 1) for ξ = 0. The kernelω appears in eq. (29) both explicitly and also within the momentum-independent expressions Next we need the entropy of the Gaussian measure, which can be calculated in the same way as for the φ 4 case (cf. appendix B). The additional colour indices on the kernelω in the measure simply yield an overall factor of N 2 − 1, so that where the trace tr is with respect to the Lorentz indices. For ξ = 0 we can now use the following The last term is ω-independent and may thus be dropped, so that 5 This is the full entropy of the measure eq. (27) with A = 0 and the symmetric kernelω. As pointed out earlier, the free action in the YM case must, however, be based on the relative 5 The formula also holds for Landau gauge ξ = 0 if the factor of d is replaced by (d − 1). entropy W eq. (12), which differs from eq. (34) by the expectation value of the Faddeev-Popov We will later treat the last term in curvature approximation, but for now we keep it general.
Eventually, we find the free action as the difference between the average action eq. (29) and the relative entropy eq. (35), B. The gap equation Our task is to minimize the free action eq. (36) with respect to the kernel ω(k). From eq. (24) and the fact that the curvature χ(k) is, in principle, a ω-dependent quantity, we have The second term describes the implicit change of the curvature with the variation kernel. While an integral type of equation for this quantity can, in principle, be written down, it represents a higher order effect and will thus be neglected within the present approximation scheme. It also vanishes for α = 1 2 . Since we will later find that the effective action depends only onω and is thus independent of α we can safely put α = 1 2 . The remaining derivative acting on eq. (36) gives the gap equation in the form For the last term, we resort again to the curvature approximation eq. (22). With the correlators from eq. (28), we have Taking the derivative w.r.t.ω(k) and neglecting again the implicit dependence of χ onω, we We can now use this result in eq. (38) and finally obtain the gap equation in the form When the solution of eq. (41) is used in the traditional way to compute Schwinger functions and eventually the effective action, the result must be because the gluon propagator at A = 0 isω −1 according to eq. (28), and the functional approach has no higher vertices.
Due to the Lorentz invariance of our approach at A = 0, the kernelω(p) is a function of |p| only and does not single out a direction, so that The gap equation (41) can now be written in a very compact form, where the gluon mass is dynamically generated through the non-linear integral equation, The gap equation (43) is very transparent: the dynamical mass M 2 is generated from the 4-gluon vertex, while the curvature χ(k) describes the coupling to the Faddeev-Popov ghost fields. As we will see shortly χ(k) is just the ghost loop, see eq. (51) below.
In fact the gap equation can be interpreted as the dispersion relation of a relativistic particle with mass M and a self-energy given by the curvature χ(k).
Let us finally emphasize that the functional approach and the gap equation (43) can be interpreted in a slightly different way: If we go back to the general definition eq. (5) of the effective action, but this time constrain the gluon propagator A a µ (−p) A b ν (p) = δ ab t µν (p)ω(p) −1 instead of the gluon field itself, we can set A = 0 in our trial measure eq. (27) and interprete the variational parameterω as the classical value for the inverse gluon propagator. The free action F (ω) from eq. (36) therefore coincides with the effective action Γ(ω) for the (inverse) gluon propagator. This entails that the optimal value forω is given by δΓ/δω = δF (ω)/δω = 0, which is exactly eq. (41) and hence the gap eq. (43). Thus, the gap equation yields the best match with the exact gluon propagator (in the sense of the effective action) which can be achieved within our variational ansatz. No such argument exists for the vertex corrections in the linear response approach, i.e. while this formulation is able to produce radiative corrections for higher-order Green functions (even with a Gaussian measure), a realistic description of higher vertices presumably requires to go beyond the Gaussian ansatz [38]. We will therefore restrict our investigations to the propagators of the theory.

C. Ghost DSE and the curvature
The gap equation (43) contains the curvature χ(k) (21), which is nothing but the ghost loop and will be calculated below. The ghost propagator is the expectation value (in our trial measure) of the inverse Faddeev-Popov operator where G 0 = − is the free ghost propgator and h = gÂ∂ describes the interaction with the gluon. 6 From the usual resolvent identities, we obtain first G = G 0 + G 0 h G and thus where we have introduced the ghost self energy −Σ = G −1 − G −1 0 . Following Ref. [23] one can derive for Σ the expression which involves the gluon propagator D from eq. (25) as well as the free (Γ 0 ) and full ghost-gluon vertex Γ defined by From Dyson-Schwinger and flow equation approaches to Yang-Mills theory in Landau gauge, it is well known that the dressing of the full ghost-gluon vertex is a subleading effect, so that the full vertex in the last equation can be replaced by the bare one (rainbow-ladder approximation).
This has the advantage that the ghost self-energy Σ (and thus the full ghost propagator G) can be expressed through the kernelω alone. To do so, we introduce the ghost form factor η via Note that both G 0 and G are colour diagonal, so that η = η(x, y) has no colour index. Using eq. (46), it is easy to see that η −1 = 1 − Σ G 0 . If we now use the rainbow-ladder approximation in eq. (47) for Σ and re-express G through the form factor via eq. (49), we obtain a closed integral equation for the ghost form factor in momentum space, where the explicit form eq. (25) of the gluon propagator was inserted.
As for the curvature eq. (21), taking two functional derivatives of (ln J ) = Tr ln G −1 and using eq. (48) yields Using the rainbow-ladder approximation again and contracting Lorentz indices and colours, we can express the scalar curvature eq. (23) in momentum space through the ghost form factor, The DSE (50) for the ghost form factor depends explicitly on the kernelω of the trial measure

A. Counterterms
To complete our analysis, we have to determine the high momentum behaviour of our Green functions and renormalise the corresponding integral equations. We begin with the ghost DSE (50). For large momenta,ω(k) ∼ k 2 and η(k) ∼ 1 (up to logarithmic corrections), so that dimensional analysis implies for the logarithmic divergence in I η (k) (50) where Λ is a suitable UV cutoff and a 0 a finite numerical constant that depends on the cutoff procedure. Since G 0 (k) = 1/k 2 , the corresponding counterterm δΣ for the ghost self-energy is proportional to k 2 , i.e. it is a ghost field renormalisation. To the given loop order, this field renormalisation leads to 1 η(k) = 1 − N g 2 I η (k) − δZ c + two loops .
In terms of explicit ghost fields, the counterterm δZ c = (δΣ)G 0 (or δΣ = k 2 δZ c ) would hence correspond to the local expression in the exponent of the trial measure (27). We can now adjust the finite pieces in δZ c such that where η(µ) is an arbitrary finite constant, because eq. (54) is independent of k and the (logarithmic) Λ-divergences on both sides agree. This prescription leads directly to the renormalised ghost equation which could also be obtained by simply subtracting the bare equation (50) at k = µ. Notice that eq. (54) may be a rather unusual field normalisation, but any other prescription for Z c can only differ by a finite constant. Notice also that the renormalised eq. (55) is independent of g, as can be seen e.g. by rescaling η → η ≡ g η.
Next, we study the mass and curvature equations (44) and (52), respectively, which are quadratically divergent by power counting, with numerical factors a i , b i and c that depend on the cutoff procedure. 7 The subtraction is a bit more complicated in this case because of the sub-leading logarithmic divergence. We begin by subtracting the k-independent contributions with counterterms for M 2 and χ(k), respectively, where χ fin is an arbitrary finite constant. The mass equation (44) is now finite with a dynamically induced mass M 2 whose absolute value is undetermined because of the finite pieces in the corresponding counterterm δM 2 1 . As a consequence, we do not need to solve eq. (44) but rather take M 2 as a finite parameter that can be chosen at will. The subtracted curvature equation takes the form This is not yet finite because the difference of the two integrals contains the subleading logarithmic divergence (cf. eq. (56)) We must therefore add a second, k-dependent counterterm δχ 2 (k) which equals the negative of the divergence on the rhs of eq. (59). In order to associate these subtractions with local terms in the exponent of the trial measure eq. (27), we isolate the k-dependent pieces in δχ 2 (k) and write the total curvature counterterm as δχ(k) ≡ δχ 1 + δχ 2 (k) = δχ 0 + k 2 δZ A with where z(µ) is again an arbitrary finite piece that can be added to the counterterm. If we now replace χ(k) → χ(k) + δχ(k) as well as M 2 → M 2 + δM 2 1 , eqs. (44) and (52) will be finite. From the gap equation (43), it is clear that the k-independent subtractions can be combined to a gluon mass counterterm δM 2 = δχ 0 + δM 2 1 , while the k-dependent counter term k 2 δZ A rescales the kinetic energy, i.e. it represents a gluon field renormalisation. These subtractions correspond to local counterterms in the exponent of the trial measure (27) (with A = 0). The renormalised curvature equation is now where χ(µ) = χ fin + z(µ) µ 2 is finite.
Eq. (62) is not very suitable for numerical investigations. To put it in a manageable form, we need to determine the constant c and devise a way to perform the necessary subtractions under the integral I χ (k). Unfortunately, I χ (k) and hence the factor c depend on the ghost profile η(k) and its derivatives, which are only known numerically. Thus, we proceed as follows: At any stage renormalisation counter term renormalisation constant gluon field δZ A z(µ) ghost field δZ c η(µ) in the iterative solution of the integral equation system, we perform the angular integration in I χ (k) with the current solution of the ghost form factor η(q) to find an expression with a complicated function f only known numerically. Our construction of the counterterms above translates into the asymptotics This can be verified numerically: In the left panel of fig. 1 as a function of q and observe that it approaches a constant value 8 for large q. The value of that constant depends on k in the expected way: In the right panel of fig. 1, we have plotted the value q as a function of k and observe that this quantity is independent of k, as expected from eq. (64).
The constant value is exactly the factor c used in eq. (62) above. Thus, at any stage of the iterative solution with the current form of η(k), we first perform the angle integrations in I χ (k) to obtain the integrand f (k, q) in eq. (63) and then evaluate the curvature χ(k) from eq. (62) written in the form This procedure determines the counterterm coefficient c from eq. (62) iteratively and renders the system of integral eqs. (43), (55) and (66) finite. It also leads to a stable solution which 8 There are still small deviations from a constant due to numerical issues that prevent us from going to very large cutoffs while preserving sufficient accuracy in the iterative solution of the integral equation system. However, extracting c as in eq. (65) yields finite iterative solutions that become cutoff-independent long before the numerics become delicate. is numerically independent of the cutoff Λ if the latter is chosen large enough. The system involves four renormalisation constants η(µ), χ(µ), z(µ) and M 2 which could even be defined at different scales µ. One of the constants is, however, redundant because finite changes in χ(µ) can be absorbed in M 2 as discussed above. In fact, only the gluon mass parameter 9 M 2 A (µ) ≡ M 2 + χ(µ) − z(µ) µ 2 will appear in the renormalized equations, and we can therefore drop either χ(µ) or M 2 in favour of the other without loss of generality. The remaining renormalisation constants η(µ) and z(µ) fix the prefactor of the kinetic energy for ghost and gluon, respectively, i.e. they fix the scale of the ghost and gluon field. (No vertex renormalisation is induced by the theory.) The three independent counter terms are also summarised in table I.

B. Numerical treatment
Let us briefly comment on the numerical treatment of the integral equation system, and in particular on the role of the various renormalization constants. First, we note that the coupling constant g can be eliminated from all equations by rescalingη = g η. If we insert the curvature equation (66) directly into the gap equation (43) and further introduce the gluon mass parameter M 2 A (µ) ≡ M 2 + χ(µ) − z(µ) µ 2 as before, we obtain the renormalized system 9 To avoid problems with possible infrared singularities, we renormalize the gluon propagator at µ > 0 so that M 2 A (µ) = M 2 + χ(µ) − z(µ) µ 2 is merely a renormalization constant without a direct interpretation as a mass. If the curvature (and thus the gluon propagator) happens to be finite at k = 0, the intercept can be interpreted as a (constitutent) gluon mass.
Recall that the coefficient c in the first equation is determined during the iterative solution of the system such that it becomes cutoff-independent.
Ideally, we would choose the renormalization point µ for the gluon and ghost field at a large Euclidean scale µ ≫ 1 far away from possible singularities, where asymptotic freedom provides natural renormalization conditions. Unfortunately, such a procedure does not lead to a stable numerical solution in the deep infrared, because the system has a family of solutions which are qualitatively different at low momenta, but cannot be discriminated in the ultraviolett. In order to stabilize the integration and bring the distinct solutions to the fore we must, at least for the ghost equation, choose a renormalization point in the deep infrared, or even at µ = 0, which we did in all numerical investigations. The gluon renormalization scale µ, by contrast, can be chosen finite and arbitrary, and we use it to rationalize all dimensionfull quantities within our numerical treatment. Thus, we find that the shape of the rescaled (dimensionless) integral equation system only depends on the three dimensionless parametersη(0) −1 , z(µ), and M A (µ)/µ.
We have already mentioned above that the first of these three parameters,η(0) −1 , discriminates between the scaling and decoupling solutions. As for the gluon mass parameter M A (µ)/µ, it has no effect on the infrared behaviour of the scaling solution, while it determines the k = 0 limit of the gluon propagator for the decoupling solution, i.e. the constitutent gluon mass. Both renormalisation parameters have negligable effect on the deep infrared behaviour but determine the momentum scale at which the transition to the infrared behaviour sets in.
Finally, the gluon field renormalization factor z(µ) changes the overall size of the kernelω(k).
More precisely, any change z(µ) → z ′ (µ) for the scaling type of solution simply leads to an overall rescalingω For the decoupling solution, this simple rescaling also holds if we simultaneously change the finite renormalization constantη −1 (0) as in eq. (68).

A. Infrared analysis
The infrared analysis of the renormalized system of integral equations (67) can be carried out much as in the Coulomb gauge case [44,45], see also Ref. [2] and [46]. If we assume a power-like behavior in the infrared,ω the non-renormalisation of the ghost-gluon vertex [41] implies the sum rule Furthermore, ifω is infrared divergent it follows from the gap equation (43) that the curvature χ(k) (52) also has the same infrared exponent α. Depending on the choice of our finite renormalization constants we obtain two different types of solutions: (i) Implementing the horizon condition [47] η −1 (0) = 0 we find the so-called critical or scaling solution with the infrared exponents which both entail α > 0, i.e. an infrared vanishing gluon propagator. These are the same infrared exponents found from the one loop DSEs in Landau gauge [2]. Only the first of these scaling solutions could be found in our numerical calculation, cf. sect. VI B. It should also be mentioned that the value of the gluon mass parameter M 2 A (µ) is irrelevant, as long as the curvature χ(k) is infrared divergent, i.e. as long as β > 1/2, which is the case for both exponents given in eq. (71). These are the so-called subcritical or decoupling solutions already found from the DSEs in Landau gauge [2]. The interceptω(0) = M 2 + χ(0) can be interpreted as a constituent gluon mass.
In the special case M 2 A (µ) = 0 there is also a solution with an infrared finite ghost form factor η(k) and an infrared vanishing gluon kernelω(k), which formally obeys the sum rule with α = −1 and β = 0. This solution is, however, definitely ruled out by all existing lattice data, in which the ghost form factor always diverges at k → 0. This leaves us with just the two type of solutions listed above. In the next subsection we compare these solution with high-precision lattice results.

B. Numerical results
As discussed above the discriminating criterion for the two type of solutions is whether we choose the ghost form factor η(0) finite (decoupling solutions) or infinite (scaling solution). Let us briefly describe our numerical findings for both type of solutions:  figure 2, the critical solution is compared to high-precision lattice data taken from Ref. [7]. It is obvious that the scaling type of solution describes the UV behaviour fairly well, but it severly deviates in the deep infrared.
Decoupling solutions: These solutions have no scaling behaviour, and instead exhibit the emergence of a soft BRST breaking mass scale in the deep infrared, so that all Green's functions remain finite at k → 0. In this case, the renormalization constant M 2 A (µ) obviously matters as it dominates the gluon propagator in the deep infrared. In addition, we can also adjust the intercept η(0) for the ghost form factor, and the overall scale z(µ) of the gluon propagator at some non-zero scale µ > 0. (The combination of M 2 A (µ) and z(µ) determines the slope of the gluon propagator in the transition region around k = µ.) The decoupling solution is shown in figure 3, along with lattice data [7] for comparison. As can be seen from the plot, this type of solution is in good agreement with the lattice data for the entire momentum range. In particular, the agreement with the ghost data is almost perfect. 10 favourably with both solutions in the UV and into the transition region, although detailed studies prefer the subcritical (decoupling) solution in the infrared. Our numerical treatment is on par with the best analytical studies, and our decoupling solution, in particular, agrees very well with the available lattice data over the entire momentum range. In addition, the variation principle used here offers a physical criterion to distinguish between the two type of solutions which is not based on arbitrary boundary conditions, namely the solution with the lowest free action must be realized. This question will be studied in a future investigation.
The method presented here works with Euclidean path integrals only, so that it naturally generalises to all extensions that can be formulated within a path integral. In particular, we can include fermions and study finite temperatures and chemical potentials without conceptual problems. These issues will also be subject to future work.
To see this, it is sufficient to find a current j = j ϕ for which the bound in eq. (A1) is saturated, because then the upper bound in eq. (A1) is a maximum (and hence a supremum). From the derivation of eq. (A1), it is clear that j ϕ should be chosen such as to obey the constraint φ µ j = ϕ for the Gibbs measure µ j that solves the minimisation problem on the rhs of eq. (A1). This is because the unconstrained minimum happens to obey the constraint and is therfore also the constrained minimum, i.e. the inequality in the third line of eq. (A1) becomes an equality. The relevant equation φ µ j = ϕ for j = j ϕ is, however, merely the extremality condition which we assume to always have a solution. 11 Thus, the effective action Γ(ϕ) defined by the variation principle (6) The corresponding value of the minimum in eq. (A1) is then given by the modified partition function with a linear term in the action, Clearly, this identifies W (j) as the usual generating functional of connected Green's functions and hence the Legendre transform eq. (A2) as the generator of 1PI proper functions. The inverse Legendre transformation implies δΓ/δϕ(x) = j ϕ (x). Finally, eqs. (A2) and (A7) entail that both functionals W (j) and Γ(ϕ) must be convex.
(An implicit regularisation is understood.) Using Wick's theorem for the correlators, all dependence on the classical field ϕ drops out, 13 The free action (2) for the Gaussian measure in φ 4 theory is now F (ω, ϕ) = S µ − W (µ) with eqs. (B4) and (B5) providing the details. At this point, the functional and linear response approach start to deviate.

The functional approach
As explained in section II, we simply take the unconstrained free action F (ω, ϕ = 0) and vary w.r.t. the kernel ω, or rather, its inverse. After a short calculation, one finds Due to the absence of a classical field, the solution of this gap equation in the functional approach has global translation and rotation invariance, i.e. the kernel has the Fourier representation where ω(k) depends only on the modulus k = |k µ |. The gap equation in momentum space takes the simple form The solution is the covariant dispersion relation for a free boson of a dynamically generated mass, where M is implicitly determined by the non-linear equation Finally, we can now compute all connected Green functions of the solving Gaussian measure, construct a generating functional W (j) and Legendre transform to obtain the effective action Γ(ϕ). Obviously, all proper vertices vanish and we get It should finally be noted that the p 0 -integration in the covariant mass equation (B10) can be performed explicitly. For d = 4, we obtain which coincides exactly with the result of the traditional variational approach in the Hamiltonian picture based on a Gaussian wave functional.

The linear response approach
In this case, we have to retain the classical field in order to comply with the constraint φ = ϕ.
From the variation of the constrained free action F (ω, ϕ), we obtain the modified gap equation As expected in the linear response approach, the optimal kernel ω ϕ determined from this equation depends implicitly on the classical field ϕ, while the effective action is simply Γ(ϕ) = F (ω ϕ , ϕ).
Since the external fields are arbitrary, the implicit ϕ-dependence of the optimal kernel ω ϕ spoils the translational or rotational invariance. As a consequence, we cannot go to momentum space as we did in eq. (B7) within the functional approach, which complicates the solution of eq. (B13) considerably.
One possible approach is to attempt a solution of eq. (B13) by expanding in powers of the classical field. Upon comparing the two gap equations (B6) and (B13), we conclude that where ω on the rhs is the translationally invariant solution of eq. (B6) found earlier. If we insert the expansion back into the free action, we find the effective action Notice that the quadratic pieces in Γ are identical for the functional and linear response approach. This is generally true because the Gaussian ansatz for our trial measure is already quadratic in the fields, and the kernels of the two approaches agree at ϕ = 0. Thus we arrive at the important conclusion: For Gaussian trial measures, the functional and linear response approach always lead to the same quadratic pieces in the effective action, i.e. to the same propagators.
The average action follows by inserting the correlators (28) into the classical action. After some lengthy calculation, we obtain S gf µ =S gf (A) + 1 2 d(x, y)ω aa µν (x, y) − x δ µν + (1 − ξ −1 ) ∂ x µ ∂ x ν δ(x, y) To find the free action, we have to subtract the relative entropy, cf. eqs. (32) and (39) It is now straightforward, though very cumbersome, to derive the gap equation by minimising F = S µ − W w.r.t. the kernelω, or rather its inverse. The result can be put in the form ω ab µν (x, y) = δ ab − x δ µν + (1 − ξ) −1 ∂ x µ ∂ x ν δ(x, y) + χ ab µν (x, y)+ Notice that the derivatives in the third line also act on the δ-function, except for the first term in the curly brackets, and the curvature also enters with its full colour and Lorentz structure.
Upon setting A = 0, only the first two lines on the rhs of eq. (C3) survive. This is just the position space equivalent to the gap equation (41), to which it reduces when the symmetric form eq. (15) is assumed forω and χ. An exact solution of eq. (C3) combined with the full curvature equation is not feasable. However, the A-dependent terms on the rhs are independent of the kernelω and thus act like an inhomogenity. This suggests and iterative solution by expanding