Dynamically massive linear covariant gauges: setup and ﬁrst results

We discuss the possibility to obtain a massive Landau gauge, based on the local composite operator (LCO) effective action framework combined with the Zimmerman reduction of couplings prescription. As a way to deal with the gauge ambiguity, we check that the ghost propagator remains positive, a necessary condition for gluon ﬁeld conﬁgurations beyond the Gribov region to be negligible. We pay attention to the BRST invariance of the construction, allowing for a future generalization to a class of massive linear covariant gauges. As a litmus test, we compare our predictions to the lattice data for the two-point functions in Landau gauge introducing the “Dynamically Infrared-Safe” renormalization scheme, including the renormalization group optimization of both the gap equation and the two-point functions. We also discuss the relation to and differences with the Curci-Ferrari model, which usefulness in providing an effective perturbative description of non-perturbative Yang-Mills theories became clear during recent years.


I. INTRODUCTION
Recently, the Curci-Ferrari model [1,2] has witnessed a revived interest in works like [3][4][5][6][7][8][9][10][11] thanks to its capability to describe rather well the n-point functions of gauge theories in the Landau gauge, next to allowing a perturbative sneak peek into the phase diagram, [12][13][14].As of now, however, the success in matching the model to lattice data relies on the fitting of both the gauge coupling g and the Curci-Ferrari mass m [4,5,8].
A possible route towards a first principle derivation of the model, including a determination of m from the sole knowledge of g at a given scale, was formulated in [15] based on a weighing over the Gribov copies.It was implemented successfully in a class of non-linear gauges that contains the Landau gauge as a limiting case [16].Unfortunately, the dynamical mass generation mechanism identified in this reference fails precisely in this limit.A more recent attempt was done in [17], directly in the Landau gauge.It was found that the system exhibits two phases, one of which corresponds to a massive implementation of the Landau gauge bearing some resemblance with the Curci-Ferrari model, with however gapped ghost degrees of freedom.As discussed in [17], the presence of massive ghosts is not incompatible with the lattice results for the latter are not a direct measurement of the ghost propagator but rather the averaging of the Faddeev-Popov operator −∂ µ D µ which remains massless in both of the above mentioned phases.A more serious difficulty is, however, that the gluon mass identified in [17] appears as a mere gauge-fixing parameter whose value is not determined in terms of g.Although it remains yet to be seen how much the correlators computed within this approach are sensitive to this gauge-fixing parameter, this could potentially compromise the comparison to lattice data.
On the other hand, it has been known for some time that dimension two condensates such as A a µ A a µ can be dynamically generated in the Landau gauge via the local composite operator (LCO) formalism [18] and the minimization of the vacuum energy.A natural question that emerges is then what connection do these condensates bear with the Curci-Ferrari model (or similar approaches), and to which extent do they allow one to reproduce the Landau gauge correlators evaluated on the lattice.The current note aims at investigating these questions.Since our goal is to eventually extend the approach to covariant gauges, we will pay particular attention to the BRST invariance of the construction.
Two dimensional condensates are of course not free of ambiguities due to their composite nature that requires extra renormalization.However, since the correlators of the primordial fields do not depend on this additional subtractions, at least not at an exact level, it is possible to exploit the reduction of coupling technique [19,20] in order to fix this arbitrariness in some way, as we recall below.In fact, one could envisage using similar ideas to fix the arbitrariness in the choice of the gluon mass in the approach of [17], at least for the evaluation of physical observables.

II. THE BRST INVARIANT CONDENSATE A. BRST invariant gauge field
Our starting point is the Euclidean Yang-Mills action in d = 4 − ε dimensions supplemented with a linear covariant gauge fixing, S with ν the non-Abelian fieldstrength tensor, and D ab µ = δ ab ∂ µ + g f acb A c µ the covariant derivative in the adjoint representation.Since our choice is to preserve BRST invariance at each step, we should only consider BRST invariant gluon mass operators, constructed out of a BRST invariant version of the gauge field.To this purpose, we insert into the corresponding path integral a unity 1 = N DξDτD η Dηe −S 1 det(Λ(ξ)), with and the appropriate normalization being collected in N .Here, we introduced the local but non-polynomial composite field with h = e igξ = e igξ a t a , where the t a denote the generators of the su(N) algebra and the ξ a are akin to Stueckelberg fields.The fields ηa and η a are additional (anti-commuting) ghosts that, together with the ξ-dependent determinant det(Λ(ξ)), account for the Jacobian arising from the change of variables ξ → A h , which is itself needed in order to treat the functional distribution δ(∂ µ A h µ ) that appears after integration over τ.As we will show in Appendix A, so long as the theory is defined in dimensional regularization, the determinant gives no contribution to the partition function at any fixed order in perturbation theory.For this reason, in what follows we will drop det(Λ(ξ)) and write the unity in the form 1 = N DξDτD η Dηe −S 1 .
We stress that, when writing equations such as (1) or (2), we are disregarding the presence of Gribov copies.The justification is two-fold.First, in this work, we restrict to perturbation theory for the evaluation of both the correlation functions and the vacuum energy.Second, we will soon check that the dynamically generated condensate is such that the ghost propagator remains positive, a necessary condition for the functional integral to be dominated by configurations within the first Gribov region.
The action (2) is here used as a way to treat, within a local setup, a BRST-invariant quantity A h µ that becomes non-local on-shell.The non-local on-shell nature of A h µ becomes explicit after one solves the condition ∂ µ A h µ = 0 iteratively for ξ using Indeed, this leads to the infinite series of terms, ∂A ∂ 2 , ∂A + . . .(7) which eventually gives the transverse on-shell expression with It can be shown that the on-shell expression (9) is gauge/BRST invariant order per order.We will have nothing to say about large gauge transformations.At the level of the off-shell/local formulation (2), the BRST invariance is explicit if one supplements the usual BRST transformation with the following transformations for the remaining fields: The BRST transformation remains nilpotent (s 2 = 0) over the extended field space and it is immediate to check that which is nothing but the infinitesimal version of For later purpose, we also mention that A h µ becomes A µ in the Landau gauge (at least perturbatively) as it is easily checked from (9) using ∂ µ A µ = 0. We refer to e.g.[21][22][23] for more details, including the connection of the representation (9) with the (local) minimization of the ℓ 2 -norm of the gauge field.

B. BRST invariant condensate
Having gone through this preparatory phase, we are now ready to investigate the vacuum structure of the local action S (2) Y M + S 1 , which is perturbatively equivalent to (1).In particular, we shall analyze the possible formation of a BRST-invariant condensate A h µ A h µ .To this purpose, we couple the corresponding gluon mass operator to a source J: 1

S
(2) Here, we introduced the relevant renormalization factors and counterterms in the J-dependent piece of the action, as these will concern us most here, 2 next to the necessary powers of µ to ensure that dim ζ = 0 (for later use, we also note that dim J = 2).More precisely, starting from bare fields, parameters and sources, which we denote by a subscript b, we write from which we deduce that the renormalization of the source is Similarly, The (pure vacuum) counterterm ∝ δζJ 2 is necessary to remove the vacuum divergences in the generating functional with W [J] = − ln D fields e −S J , hence its appearance as an additive renormalization.Let us stress that these vacuum divergences do not affect the divergences appearing in correlation functions with at most one insertion of A h A h .Therefore, we can (and will) consider only renormalization schemes where all renormalization factors, including Z 2 , are ζ-independent.
In general, from ( 16) it follows that 1 The method presented here was first worked out in [24], see also [18,25] and [26].We will however slightly adapt the discussion to point out some new, yet unnoticed, features, which also facilitate the interpretation. 2The other corresponding factors in S where µ ∂J ∂µ = γ J J. We immediately discarded terms that vanish in the ε → 0 + limit.We will recall below that, in principle, it is possible to choose δζ merely proportional to ζ such that ζ + δζ = Z ζ ζ, henceforth implementing multiplicative renormalization also for the vacuum divergences.
We note that from (17), or directly from (16), it follows that a finite shift implies that making ζ + δζ invariant.The shift δζ f in can depend on all other variables.The invariance of ζ + δζ under such shifts was first discussed in [25].
From the generating functional W [J], one can introduce the field which becomes the argument of the effective action Γ(σ) ≡ W (J) − d d xJσ, with δΓ/δσ = −J.As usual, the benefit of such a construct is that it allows one to access the zero source limit of σ, and therefore A h µ A h µ J→0 , from the minimization of Γ[σ].In particular, a non-trivial minimum is tantamount to the dynamical generation of a condensate A h µ A h µ J→0 = 0. We mention here that the non-positivity of the integration measure associated to the action (1) could potentially jeopardize the usual relation between the limit of zero sources and the minimization of Γ[σ].However, we shall later verify that the dynamically generated condensate is such that the ghost propagator remains positive, suggesting that field configurations with negative measure have a subdominant effect and, therefore, that the minimization prescription can be used.Strictly speaking, a similar check should be done with the new ghost fields η and η.
To actually compute the effective action, it is computationally simplest to rely on Jackiw's background field method [27,28].In the present context, this is not easily done a priori due to the coupling of J to a composite operator and to the presence of the quadratic term ∝ J 2 .We can however easily remedy this situation following [18,24].We insert a unity 1 = N Dσ e −S σ , with The last two terms of the third line cancel exactly the two J-dependent terms in (13) hereby defining a new sourced action

S
(3) (3) Y M + S 1 + S 2 , and The source now couples linearly to a primary field σ (whose expectation value is of course (19)) and the background field method can be implemented as usual.Of course, integrating exactly over σ, working with S Y M will give completely equivalent result as with S (2) Y M .The situation will only get interesting if the dynamics of the theory would prefer to assign a non-vanishing vacuum expectation value to σ.This is a possibility that we now investigate.Before doing so, we also notice that the bare action only depends on the combination 3 ζ + δζ, which we already argued to be independent of the finite parts in δζ.Without loss of generality, we can thus renormalize the vacuum divergences in a computationally efficient scheme as MS.
To do so, we first notice that, given the BRST invariance of both the action and the mass operator A h µ A h µ as well as the fact that the α-dependent part of the action is BRST-exact in the limit of zero sources, the expectation value of σ does not depend on α in this limit. 4Therefore, we can (and will) choose to work in the Landau gauge, α → 0, in which case the (τ, ξ)-integration can also be done exactly, leading to the on-shell identification A h → A [29].At one-loop order, one obtains where, as usual, we have divided the effective action by the space-time volume to compute the effective potential V (σ) and we have included an extra factor µ ε to ensure that dimV = 4.We have also treated δZ 2 ≡ Z 2 − 1 and δZ ζ ≡ Z ζ − 1 = δζ/ζ as higher loop corrections, neglecting them in the one-loop term, and expanding to first order in δζ in the tree-level term.After explicitly computing the integral and absorbing the divergence in δζ using the minimal subtraction scheme, one arrives at the expression where we have defined

C. Reduction of couplings
Various remarks are in order at this point.First, the trivial solution m 2 = 0 is always a maximum of the potential since V ′′ (m 2 → 0) = −∞.This shows that, in the present approach and at the present order of evaluation, a condensate A h µ A h µ J→0 is dynamically generated (independently of the value of ζ).
Beyond this proof of existence, the next pressing question is the size of the condensate.Here, however, we face a serious problem: while the running of ζ with μ is entirely fixed from the renormalization factor, its value at a chosen initial scale μ0 is arbitrary and impacts directly on the size of the condensate.On the other hand, it is easily shown that δW /δζ ∝ J 2 .This means that, were we to work exactly, ζ should not influence any quantity with less than two J-derivatives, in the limit of zero sources, or, in other words, at the minimum of the effective action.This includes the correlation functions for the primary fields but also the condensate itself A h µ A h µ J→0 .Of course, at a given order of approximation, a spurious dependence with respect to ζ is to be expected but the previous argument shows that, given a prescription for choosing ζ, we can test a priori how the ζ-independence (re)emerges as the approximation is improved.Let us now discuss two possible prescriptions that could be used.
A first possibility would be to impose the ζ-independence of certain quantities such as the condensate or the value of the potential at the minimum.At the present order, this is not very useful because 1) the constraint does not fix any particular value of ζ and 2) it leads to m 2 = 0 which we have already argued to correspond to an unstable state.It would be interesting to see how this constraint is modified at the next order of approximation.We leave this interesting question for a future investigation.
Another possibility is to follow the prescription of [18,24].Indeed, despite the presence of two coupling constants, g 2 and ζ, g 2 runs as usual, that is separately from ζ.Moreover, at an exact level, the value of ζ does not affect the quantities that we are after, namely the correlation functions for the primary fields.This makes our situation suitable for the Zimmermann reduction of couplings program, [19], see also [20] for a recent overview, in which case one coupling (here ζ) is re-expressed as a series in the other (here g 2 ), so that the running of ζ controlled by γ ζ (g 2 ) is then automatically satisfied.This selects one consistent coupling ζ(g 2 ) from a whole space of couplings ζ, and it is also the one (unique) choice compatible with multiplicative renormalizability, ζ + δζ = Z ζ δζ.This approach was also applied to the Gross-Neveu model in [30], reporting very good agreement with the exactly known mass gap in this toy model.
In the MS scheme, one finds [3,18], which is a particular solution of which evidently draws from (17).
From our point of view, this choice is consistent with the perturbative approach followed in this work, where any quantity is assumed to admit an expansion (be it a Laurent expansion) in powers of g 2 .

III. TWO-POINT FUNCTIONS
Let us now study how the dynamically generated condensate affects the Landau gauge two-point correlation functions.Rewriting the field σ in (22) as it vacuum expectation value which we also denote σ and a fluctuating part δσ, we find that the free gluon propagator is similar to the one in the Curci-Ferrari model The ghost propagator is simply G(p) = δ ab /p 2 , while the δσ propagator is ζ/µ ε .Moreover, in addition to the usual Landau gauge vertices, we have a AAδσ vertex and a new AAAA vertex − µ ε 4!ζ δ ab δ cd δ µν δ ρσ + δ ac δ bd δ µρ δ νσ + δ ad δ bc δ µσ δ νρ .
(32) We mention that, strictly speaking, ζ −1 , and therefore m 2 , have a perturbative expansion in terms of ζ 0 , ζ 1 , . . .which one needs to take into account in order to evaluate the correlation functions at a given order.One could define the Feynman rules in terms of this expanded parameters.However, it is more convenient, and equivalent, to consider the Feynman rules in terms of ζ −1 and m 2 and re-expand them only at the end of the calculation, when necessary.

A. Ghost propagator and Gribov horizon
Using the above derived Feynman rules, we find that the one-loop ghost propagator coincides with the one computed in the CF model.We can use this remark to elucidate the role of the dynamical condensate within the context of Gribov's construction [32] to avoid multiple solutions to the Landau gauge condition.Following Gribov's original setup [32], the ambiguity related to the presence of Gribov copies is handled by restricting the domain of integration in the functional integral to the Gribov region As it is apparent from the definition of the region Ω, the ghost propagator ca c b p , i.e. the inverse of M ab , remains positive within Ω.The positiveness of the ghost propagator is thus a necessary condition for any approach to be compatible with a restriction to the Gribov region Ω, and we can check to what extent our gluon mass scale is consistent with such condition.
We parameterize the ghost two-point vertex as where is the ghost self-energy and we have factored out the trivial color structure δ ab .At one-loop order, one finds -modulo renormalization- where, to the present accuracy, m 2 can be expanded to leading order in g 2 that is replaced by m 2 0 (which we keep denoting m 2 for simplicity from here onwards).As σ gh (p 2 ) is a decreasing function of the momentum square p 2 [33], the positivity of the ghost propagator will be ensured by demanding that σ gh (0) < 1.Using dimensional regularization in the MS scheme, 5 one finds where we have reparameterized the coupling by defining This leads to in the zero-momentum limit.The positivity condition translates then into On the other hand, the minimum of the potential ( 24) is located at It is easily checked that this value obeys the positivity bound (38) for any value of λ, clarifying then the role the dynamical condensate plays in relationship to the issue of gauge copies.

B. Gluon propagator
Similarly, the gluon two-point vertex reads

Γ
(2) where we have again factored out the trivial color structure δ ab and it is understood that m 2 and ζ need to be expanded to the appropriate order, namely to leading order, except for the tree-level mass term m 2 that needs to be expanded to nextto-leading order according to (28).Up to the combination Z 2 − Z ζ whose relation to the mass counterterm in the CF model we have not worked out yet, the first line is the oneloop Γ (2) AA as computed in the CF model with mass m 2 .The two other contributions correspond to two new diagrams that arise from the vertices (31) and (32).Now, the terms involving D µν cancel between the last two lines of ( 40) and we are left with Writing (23) in terms of m 2 and setting V ′ (m 2 ) = 0, we find This implies that For the above equation to make sense, we see that m 2 δZ 2 should correspond to the mass counterterm in the CF model We have checked that this is indeed the case using the value for δZ 2 given in [18].This is of course no surprise since a constant source J plays exactly the same role as the Curci-Ferrari mass and therefore Z J = Z CF m 2 , or, owing to (15), Z 2 = Z A Z CF m 2 .Interestingly, the constant ζ disappears from the gluon self-energy, see (43), once the latter is computed on the shell of the gap equation.This could have been foreseen, given that such a parameter was not present in the Faddeev-Popov action in the first place, see our earlier made comments.Nonetheless, an implicit dependence on ζ still survives via the solutions of the mass gap equation.
From (43), it is evident that, although the condensate acts as a mass for the tree-level propagator, it disappears from the tree-level term at one-loop order. 6Thus the LCO approach does not lead to the CF model but rather to something more akin to the screened massive approach of [35][36][37] where the mass disappears from the loops as the order of approximation is increased.Of course, such a reorganization of the 6 Strictly speaking, the gluon propagator is obtained after considering δσA and Γ Fortunately, using that ζ −1 ∼ g 2 , it is easily shown, however, that Γ (2) δσδσ ∼ g 2 which means that the extra term can be neglected at the present order of accuracy.We also note that Γ Aδσ vanishes in the zero-momentum limit from symmetry considerations.perturbative expansion does not necessarily lead to a trivial reformulation of the theory without condensates.In this respect, it would be interesting to push the present approach to the next order and identify quantities that remain constant and different from their corresponding predictions in the absence of condensate.We leave this interesting question for future work.Nonetheless, in Appendix B we already offer a generic argument as to why the tree level mass indeed should cancel upon using the gap equation.

C. The Dynamically Infrared-Safe (DIS) scheme and the Renormalization Group
We would now like to compare our results to lattice data.To this purpose, we need to choose a renormalization scheme that is free of Landau poles and that yields a gluon propagator which displays the desired infrared behavior.The renormalization factors we need to fix are Z A , Z c , Z λ and Z 2 .We will also redefine Z ζ to better match our scheme, rather than using MS like we did in Sec.IC.Since Z ζ ultimately only affects the solution to the gap equation, we will postpone this to the next section.
To begin with, since we are in the Landau gauge, in order to fix the coupling renormalization Z λ we can impose the Taylor condition In the Taylor scheme, the running of λ is completely determined by the anomalous dimensions γ A and γ c of the gluon and ghost propagators, where Z A and Z c themselves can be defined in the momentum subtraction (MOM) scheme, as is appropriate for a comparison with the lattice data.Namely, if we renormalize the gluon and ghost propagators at the scale µ so that (here Γ (2)⊥ AA denotes the transverse component of the gluon two-point vertex), then, going back to ( 43) and (33), where Π m 2 ⊥ CF,1ℓ is the transverse component of the CF gluon polarization [5], Observe that the first of the conditions in (47) is especially suitable to our model, which does not have a tree-level mass term in the inverse propagator.Indeed, due to the cancellations which occur in (41) as soon as the gap equation is enforced, ( 47) is favored over conditions such as Γ (2)⊥ AA (µ 2 ) = µ 2 + m 2 , which is often used for renormalizing the gluon propagator in the CF model [5].
At this stage, δZ 2 in ( 48) remains yet undetermined.Before discussing its definition, let us explore how Z 2 is relevant to renormalization, starting from the gluon mass m 2 .In the present model, we do not really have a mass renormalization factor since m 2 = −g 2 σ/ζ 0 is defined only at the renormalized level.This is at variance with the CF model, where it is found [4,5,26] that β m 2 = m 2 (γ A + γ c ).However, we can rewrite m 2 in terms of µ-independent bare quantities modulo µ-dependent renormalization factors.Explicitly, we have It follows that the µ-dependence of m 2 is encoded in the prefactor where we have used the Taylor condition, (44).In particular, the µ-dependence of m 2 is entirely governed by Z 2 and Z c , with the square mass beta function β m 2 reading Furthermore, as we have seen in Sec.IIB, Z 2 plays the same role as Z A Z CF m 2 in the CF model.Since the latter enjoys a non-nilpotent BRST symmetry which implies that Z CF m 2 Z A Z c is UV finite, we deduce that, in the present model, Z 2 Z c should also be UV finite. 7The most obvious condition to impose for defining Z 2 would then be Z 2 Z c = 1 [26].However, there are good reasons to use a slightly different definition.
To see why an alternative condition to Z 2 Z c = 1 is appropriate in the present renormalization scheme, it is instructive to compute the low-energy limit of the beta function β λ while keeping δZ 2 arbitrary.As µ 2 → 0, we have-see ( 35), ( 37) and ( 50)- where the O(µ 2 ) term in Π m 2 ⊥ CF,1ℓ contains a logarithm of µ 2 .Plugging these into ( 48) and ( 49), and using (46) while neglecting any higher-order term in the renormalization scale and in the coupling, we find that the gluon and ghost anomalous dimensions have the following asymptotic behavior: By (45), these yield For most values of δZ 2 , subject to the sole condition that γ A (or, equivalently, β λ ) do not contain divergences, we find that which is equivalent to8 This behavior is consistent with the one found for the Taylor coupling on the lattice [38].On the other hand, suppose that δZ 2 = −δZ c , as would be the case to one loop if we decided to choose our renormalization condition according to Z 2 Z c = 1.Then we would have, by ( 49) and ( 55), yielding β λ → constant × λ 2 in the low-energy limit.This is equivalent to λ(µ) ∝ 1/ ln µ 2 → 0 as µ → 0.
While having a running coupling that vanishes logarithmically at zero momentum might not be an issue in and of itself, if we insist that Z 2 Z c = 1 in the Taylor scheme, then we run into the conclusion that the gluon propagator remains massless at low energy, when improved by the methods of the Renormalization Group (RG).In order to see this, let us first write down the explicit expression for the RG-improved gluon two-point vertex in the present scheme.Denoting by Γ (2)⊥ AA (p; µ) the two-point vertex renormalized by minimal subtraction at the scale µ, one has As already observed in [39] based on the relation , The gluon anomalous dimension γ A can then be re-written as where in deriving the above equality we used the Taylor condition.We then compute Our final expression for the RG-improved gluon two-point vertex therefore reads We remark that this expression is valid in any scheme in which γ 2 = −γ c holds together with the Taylor condition.Now, since according to (57) (and by dimensional counting) using β m 2 = γ c m 2 , we find that Provided that the running coupling vanishes at zero momentum (whether logarithmically or quadratically is irrelevant), the above equation implies that m 2 (µ) saturates to a constant in the low-energy limit.Therefore, as p → 0, we find that We now see why having a logarithmically vanishing coupling spoils the infrared behavior of the RG-improved gluon propagator: if λ(p) ∼ 1/ ln p 2 , then Γ (2)⊥ AA (p) ∼ p 2 ln p 2 , implying that the propagator diverges at p = 0.
Ultimately, the reason why we found this divergence is that the specific choice Z 2 Z c = 1, which to one loop is equivalent to δZ 2 = −δZ c , kills the deep-infrared m 2 /µ 2 term in the beta function β λ .Using such a condition, however, is not at all forced on us by any consistency requirement: so long as the divergences in δZ 2 and δZ c are taken to be equal with an opposite sign, we are free to choose the finite part of Z 2 according to the needs of our scheme.Therefore, instead of picking the single value that yields a diverging propagator, we decide to use a slightly different renormalization condition for δZ 2 , namely where [δZ c ] fin. is the finite part of δZ c .In other words, we start from δZ c and subtract from it its zero-energy finite term.Note that in the MOM scheme any µ-dependent term must actually be regarded as part of the divergence, since µ is an arbitrary scale introduced by dimensional regularization which needs to be absorbed by renormalization.
Going back to (61), we see that, to leading order, the above condition is equivalent to Immediately, we find that in the scheme defined by (70), which henceforth we will refer to as the Dinamically Infrared-Safe (DIS) scheme, the beta function β λ regains its infrared m 2 /µ 2 term: As discussed earlier, this implies that λ vanishes quadratically as µ → 0. Furthermore, since the difference between δZ 2 and −δZ c is a µ-independent constant, the relation γ 2 = −γ c still holds despite being Z 2 Z c = 1.Therefore, the expression in (66) for the RG-improved gluon two-point vertex is still valid in the DIS scheme, as well as the asymptotic behavior described by (69).However, this time λ(p) ∼ p 2 , so that and the propagator saturates to a constant at zero momentum, as it should be.
In the present scheme, the gluon and ghost anomalous dimensions γ A and γ c read where t = µ 2 /m 2 .In the UV, the ordinary Yang-Mills behavior, namely is recovered.As we will see, the RG-improved gluon and ghost propagators computed from these anomalous dimen-sions are in very good agreement with the lattice data over a wide range of momenta.
Finally, let us write down an expression analogous to (66) for the RG-improved ghost two-point vertex.In general, the two-point vertex Γ (2) cc (p; µ) renormalized by minimal sub-traction at the scale µ can be written as Since in our scheme γ c = β m 2 /m 2 , we compute In particular, the ghost two-point vertex can be expressed in terms of the running gluon mass squared m 2 (p) as The above equation tells us that the ghost propagator diverges like 1/p 2 at zero momentum, and that the ghost form factor-as a function of momentum-has the same behavior as the running gluon mass squared, being essentially equal to it modulo a constant factor of m −2 (µ).
In the next section, we will redefine the gap equation in order to derive suitable initial conditions for the RG flow of the theory.These will be used to compute the RG-improved gluon and ghost propagators, which we will then compare with the lattice data in a later section.

D. Gap equation revisited and the RG flow
The gap equation allows us to fix the initial value of the gluon mass m 2 in the RG flow starting from the initial renormalization scale µ and coupling λ.Analogously as for the propagators, we will benefit from the RG invariance of the effective potential (29) to optimize it.We can first work in the MS scheme to do so.As per construction µ d dµ V = 0, or, we can use this RG equation to resum all leading logs.Therefore, we set , we get from (80) at leading order This can be rephrased as and the optimized potential becomes The above RG methodology was borrowed from [40], see also [41] for further interesting comments about RG log resummations.In Appendix C, we have taken a closer look at the effective potential.
It can be easily seen that at this order, the last factor amounts to replacing m 4 (µ) 2g 2 (µ) → m 4 (m) 2g 2 (m) , clearly showing potentially large logs are resummed, and thereby also establishing the explicit µ-independence of the improved effective potential, see also [41].
Expressing everything in terms of λ, we obtain as solution of ∂V LL ∂m = 0 in the MS scheme m sol.,MS (µ) = µe as , see e.g.[3] and (63).For later usage in Appendix C, we already mention 24(16π 2 ) 2 .In order to make use of this solution in the DIS scheme, we first need to perform a scheme conversion from MS to DIS, using the appropriate renormalization factors.Going back to (44) and (52) we see that, at a fixed renormalization scale, or, to one loop, since Z 2,MS Z c,MS = 1, Of course, to lowest order, the mass and coupling on which the renormalization factors depend can be computed in any of the two schemes.
In what follows, we will fix the renormalization scale to be equal to µ 0 = 1 GeV-i.e., the scale at which we start the renormalization group flow-and use λ DIS (µ 0 ) as our independent variable.Then we will compute λ MS (µ 0 ) as a function of λ DIS (µ 0 ) using (90), plug the result into (86) to obtain the solution of the MS resummed gap equation, and finally convert the solution to the DIS scheme using (89).Doing so yields the DIS mass m 2 DIS (µ 0 ) corresponding to the DIS coupling λ DIS (µ 0 ) by virtue of the MS gap equation.In Fig. 1 we show the solutions of the gap equation, both in the MS and in the DIS scheme, at the scale µ 0 = 1 GeV.As we can see, the conversion modifies only slightly the relation between the mass and coupling, at least for not so large values of λ.
Dropping the subscript for the DIS quantities, with the DIS m 0 = m(µ 0 ) = m sol.computed at the initial scale µ 0 as a function of the DIS coupling λ 0 = λ(µ 0 ) as detailed above, in Fig. 2 we display some of the solutions to the beta-function equations in the DIS scheme.As we can see, the running coupling has no infrared Landau pole 9 and vanishes like p 2 as p → 0, as we anticipated in the last section.The absence of Landau poles in the running coupling is a necessary condition for the self-consistency of any perturbative approach to Quantum Chromodynamics, and one which is shared by most of the models which treat the gluons as massive-see e.g.[10,37].Indeed, it is precisely the existence of a gluon mass scale that makes it possible for the running coupling 9 We checked that this behavior persists to very large values of the coupling λ 0 at the initial scale µ 0 = 1 GeV.Indeed, we could not find any value of λ 0 for which the running coupling diverges.to change its behavior and decrease as the momentum decreases.At large energies, since β λ → − 22λ 2 3 as µ → ∞ by (74) and (75), the coupling runs just like in ordinary Yang-Mills theory.
Finally, at zero momentum the running gluon mass saturates to a constant m(0), again in agreement with our analysis of Sec.IIC.Since β m 2 < 0, with m 0 = m sol.an increasing function of λ 0 , m(0) increases with the initial value of the coupling.It attains the expected order of magnitude when λ 0 0.2 − 0.3, which corresponds to values of the coupling λ 0.4 − 0.5 at the peak (equivalently, α s 1.7 − 2.1).In the UV, we find that yielding a mass that decreases like λ 9 44 (µ) ∼ [ln µ 2 ] − 9 44 at large energies, thereby restoring the massless, asymptotically free, UV limit.FIG.3: RG-improved gluon propagator (top), gluon form factor (middle) and ghost form factor (bottom) in the DIS scheme, renormalized at µ 0 = 1 GeV, together with the lattice data of [38] and analogous CF model results (red curves) for comparison, see the text for details.The initial value of the coupling λ 0 = 0.473 was obtained by fitting the combined lattice data for the gluon and ghost form factors.The initial value of the gluon mass m 0 = 0.655 GeV was computed by using the gap equation as described in Sec.IIID.

E. Comparison with the lattice data
We are now in a position to compare the predictions of our model with the lattice data, more precisely the 80 4 , β = 6.0 gluon and ghost data sets of [38], see also [42].In Fig. 3 we show the RG-improved gluon propagator, gluon form factor corresponding to the fit in Fig. 3. and ghost form factor renormalized at the scale µ 0 = 1 GeV.In order to obtain the two-point functions, we fitted the initial value of the coupling10 by the least-squares method using the combined lattice data of [38] for the gluon and ghost form factors, with the gluon mass m 0 = m sol.computed from the gap equation as a dependent input, as described in the previous section.The fit yielded a value of λ 0 = 0.473, which corresponds to α s (µ 0 ) = 1.981, m 0 = 0.655 GeV, and λ MS (µ 0 ) = 0.316.
Both the form factors turn out to be in very good agreement with the lattice data at moderate to high energiesup to 8 GeV-, despite the gluon falling slightly below the lattice in the UV, and slightly above it at intermediate energies.On the other hand, at momenta p 0.5 GeV, the RG-improved functions are suppressed with respect to their lattice analogues.This behavior is not unseen at one loop in massive expansions of Yang-Mills theory [43], at least as far as the gluon is concerned.Indeed, it is shown by the CF model itself, when the latter is fitted to the lattice data by the same procedure used for our fit 11 .Previous studies of the CF FIG.5: Running coupling corresponding to the fit in Fig. 3, rescaled by a factor of 1/2.2.Lattice data from [38].See the text for details.
model [8] suggest that this behavior could improve by going to two loops.
The running coupling λ(p) and mass m(p) given the initial value λ 0 = 0.473 and the gap equation are shown in Fig. 4. The maximum of the coupling occurs at p ≈ 0.64 GeV, where λ ≈ 0.66 (i.e.α s ≈ 2.77).At zero momentum, the running mass saturates to m(0) ≈ 0.78 GeV.
In Fig. 5 we compare the lattice data for the Taylor coupling to our running coupling α s (p).Since the fitted λ 0 (equivalently, α s (µ 0 )) is a one-loop estimate of the coupling, a rescaling of α s (p) is needed in order to match the lattice.In our case, we had to divide α s (p) by 2.2 in order to align the former both to the lattice UV tail and to the value of α s at the initial renormalization scale µ 0 .At µ 0 = 1 GeV, the rescaled value of the coupling α s is found to be 0.90.

F. Testing the stability of the DIS scheme
In order to test the stability of the DIS scheme, it is useful to extend the results presented in Secs.IIIC-IIIE to different initial renormalization scales µ 0 and to alternative renormalization parameters.
For the purposes of this section, we shall denote with ∆Z studies on the subject, but rather it was made anew starting from the combined gluon and ghost lattice form factors of [38] for the purpose of comparison.In more detail, the RG-improved CF functions were computed in the so-called infrared-safe scheme [5], rescaled so that Γ (2)⊥ AA (p)/p 2 = 1 at p = 1 GeV, like in our DIS scheme.The rescaling factor could as well be determined by leaving it as a free parameter of the fit, in which case the CF propagators would significantly improve-especially in the UV.However, we chose not to do so, in order to make the comparison to our results more immediate.the quantity defined as In our previous analysis, we took ∆Z = 5λ 8 , the latter being equal to the one-loop zero-momentum limit of the finite part of δZ c .While natural in many respects, this choice is far from unique: in Sec.IIIC we saw that any non-zero ∆Z of the form const.× λ yields a gluon propagator which saturates to a finite value as p → 0. It thus makes sense to ask whether choosing a constant other than 5/8 in ∆Z would substantially alter our results.
In Figs. 6 and 7 we answer this question by fitting the lattice data while setting ∆Z = 5λ 8 +20% −20% -that is, ∆Z = λ 2 , 3λ 4 .Additionally, we integrate the RG flow starting from three different initial renormalization scales -µ 0 = 1 GeV, 2 GeV and 5 GeV -in order to test the scale-dependence of the scheme.The fit parameters -reported in Tab.I -were obtained by adapting the procedure laid out in Sec.IIID-IIIE to the new renormalization parameters.
It is clear from the fits that increasing or decreasing the value of the coefficient of λ in ∆Z does not have a significant impact on the qualitative behavior of the scheme.Quantitatively, a change in the initial values of the coupling constant λ 0 largely compensates for the difference in the coefficients, although some deviation from our previous results can be observed in the intermediate-to high-energy regime of the gluon sector and in the low-energy regime of the ghost sector at µ 0 = 1 GeV and -less prominently -at µ 0 = 2 GeV.
As for the change of the initial renormalization scale, increasing µ 0 yields a worse match with the lattice data, especially at intermediate energies in the gluon sector and at low energies in the ghost sector.Nonetheless, we assess that the agreement between the dynamical model/DIS scheme and the lattice can be still considered satisfactory, given that the present results are obtained by a one-loop calculation and by fitting a single free parameter.Bottom: µ 0 = 5 GeV.See the text for details.

G. SU(2)
To end our overview of the dynamical model, we should mention that some preliminary tests were performed in order to evaluate whether the DIS scheme is also capable of capturing, qualitatively and quantitatively speaking, the SU(2) pure Yang-Mills dynamics.In general, as displayed in Figs. 8  and 9, these tests showed a worse agreement with the lattice data of [34,[44][45][46] in comparison to SU (3).We believe this could be mostly due to a failure of the simple one-loop approximation in the ghost sector, exemplified by the fact that we were not able to obtain a fit of the lattice data for the FIG.7: RG-improved ghost form factor in the DIS scheme at different initial renormalization scales µ 0 and for different ∆Z's.Plots as in Fig. 6.
ghost form factor alone.
In more detail, we performed simple two-parameter fits of the SU(2) gluon and ghost form factors, both separately and simultaneously, at the initial renormalization scale µ 0 = 1 GeV.The simultaneous fit of the two form factors -shown in Fig. 8 -yielded an unsatisfactory agreement with the lattice data of [34,[44][45][46], especially in the gluon sector.A good match of the gluon form factor/propagator, on the other hand, was obtained at the expense of the ghost form factor by fitting the former on its own -see Fig. 9. Interestingly, in this second case, the fitted value of the coupling constant λ(µ 0 ) = 0.35 was found to be nearly half that obtained by simultaneously fitting the two form factors -i.e., FIG.8: RG-improved SU(2) gluon propagator (top), gluon form factor (middle) and ghost form factor (bottom) in the DIS scheme, renormalized at µ 0 = 1 GeV, together with the lattice data of [34,[44][45][46] and analogous CF model results (red curves) for comparison.The gluon and the ghost form factors were fitted simultaneously, see the text for details.
λ(µ 0 ) = 0.59.As for the fit of the ghost form factor alone, our employed algorithm was not able to reach convergence and provide us with meaningful values of the parameters.
Some insight into these results can be gained by taking the Curci-Ferrari model as a reference.A reiteration of the SU(2) Curci-Ferrari fits [8] using the procedure described in Note 11 -see Figs. 8 and 9 -displays most of the features we just reported for the DIS scheme: fitting the gluon form factor alone does not yield a good agreement with the lattice data for the ghost sector, which only improves when a si- FIG.9: RG-improved SU(2) gluon propagator (top), gluon form factor (middle) and ghost form factor (bottom) in the DIS scheme, renormalized at µ 0 = 1 GeV, together with the lattice data of [34,[44][45][46] and analogous CF model results (red curves) for comparison.Fit of the gluon form factor alone, see the text for details.
multaneous fit of both the form factors is performed; this in turn requires a more-than-doubled value of the coupling constant 12 , and leads to a worse match in the gluon sector.We should note that, as far as the dynamical model/DIS scheme is concerned, this kind of behavior is unseen in SU (3), where fitting the lattice data to the gluon and/or the ghost form factors, either simultaneously or separately, yields essentially the same value of the coupling constant 13 .Despite these similarities, however, there is one major difference between the two models: just like in SU(3), instead of quickly saturating and deviating from the lattice data as in the DIS scheme, the one-loop RG-improved Curci-Ferrari ghost form factor manages to reproduce the data well down into the deep IR.
In fact, at variance with the DIS scheme, fitting the ghost form factor alone is actually possible within the CF modelat the price of a further increase in the value of the coupling constant 14 .
We can then conclude that two main factors are at play in making the DIS scheme perform worse in SU(2) than in SU(3): 1. at one loop, a sub-optimal agreement within the ghost sector makes it harder to obtain a good overall fit of the lattice data; 2. this mismatch pushes towards larger values of the coupling constant, were the one-loop approximation is less trustworthy.Since finding a good two-parameter fit of the SU(2) data was not possible within the present approach, we did not try implementing the gap equation and attempt a one-parameter fit.We leave a more in-depth analysis of the SU(2) case to a future study.

IV. CONCLUSIONS
We showed that the non-vanishing of a non-local BRST invariant mass dimension two condensate in pure Yang-Mills theory can be probed in any linear covariant gauge by minimizing an effective potential for the condensate, leading to a dynamical gluon mass which value is fixed in terms of a gap equation.For renormalization purposes of the potential, we had to introduce a novel coupling, but by resorting to a procedure known as the reduction of couplings, this parameter could be expressed as a power series in the usual coupling constant.
The computation of the potential and condensate was carried out in Landau gauge, in which case the non-local condensate reduces to just A 2 and computations drastically simplify.
The renormalization group improvement of the gluon and ghost propagators was performed in a newly introduced renormalization scheme termed the Dynamically-Infrared-Safe (DIS) scheme.In the infrared, the dynamically massive model reproduces the expected, non-perturbative behavior of pure Yang-Mills theory not only qualitatively-by the saturation of the gluon propagator at zero momentum-but also quantitatively, as demonstrated by a comparison with the lat-tice data.In the UV, where the effects of the gluon condensate are negligible as shown explicitly by renormalization group arguments, it reduces to ordinary perturbation theory.
Our results thus indicate that the BRST invariant gluon condensate is a good candidate for explaining by which mechanism dynamical mass generation occurs in the gluon sector of pure Yang-Mills theory 15 .Albeit that the eventually dynamically massive model shares many similarities with the Curci-Ferrari model, it is different as highlighted in the text.
A logical next step would be to extend the propagator and renormalization group analysis to other linear covariant gauges, in which case the Nielsen identities [53][54][55][56] might prove valuable to get a grip on the gauge (in)dependent contributions.Going beyond the Landau gauge will also require to take into proper account the then explicit non-local nature of the condensate, based on [57].
The latter is obtained by a change of variables F → ξ in the functional integral followed by the factorization of the functional determinant and the re-writing of det −∂ • D(A h ) in terms of a functional integral over a pair of ghost fields (η, η), and of δ(∂ • A h ) in terms of its τ-Fourier transform.This is analogous to the Faddeev-Popov procedure which is routinely carried out to derive the partition function of pure Yang-Mills theory in the Landau gauge, with the ghosts η and η in place of the standard ghost fields c and c, and the Fourier field τ in place of the Nakanishi-Lautrup field b.
In the context of gauge fixing, thanks to the gaugeinvariance of the partition function, one usually exchanges A h µ for A µ in (A1)-(A2) by a change of the gluon field integration variables, after which the only ξ-dependent term left in the unity is the determinant det(Λ(ξ)).It is easy to see, then, that the ξ-integral decouples from the rest of the partition function, so that the former can be absorbed into the normalization factor N .On the other hand, our introduction of the unity in Sec.IA is carried out at a stage in which the partition function is already gauge-fixed.Therefore, no change of variable A h µ → A µ can be performed, and, in a generic linear covariant gauge, the non-decoupled field ξ must be treated on the same footing as the other dynamical fields of the theory.
Nevertheless, it can still be shown that the determinant det(Λ(ξ)) does not perturbatively contribute to the n-point functions of the theory, as long as it is defined in dimensional regularization.As a consequence, when doing calculations in perturbation theory using dimensional regularization, the determinant can be suppressed by setting det(Λ(ξ)) = 1.
In order to prove our statement, we first rewrite the determinant in terms of a functional integral over a new pair of ghost fields (λ, λ), We may reexpress (A3) as det(Λ(ξ)) = DλDλ e −(I 0 +I 1 ) , where having defined The action term I 1 contains the interaction between (λ, λ) and ξ.The latter is quadratic in the ghost fields, with its ξ-dependence encoded in the function Ω ab (ξ).I 0 , on the other hand, contains the zero-order ghost propagator, which is easily seen to be Q ab (p) = δ ab in momentum space, or Q ab (x) = δ ab δ(x) in coordinate space.Now, consider the vacuum expectation value O of an operator O which does not depend on the newly-introduced fields (λ, λ).This can be computed as where the subscript 0 denotes that the average is to be taken with respect to the action I 0 plus any other (λ, λ)independent term originally present in the full action of the theory.OI n 1 0,conn.explicitly reads , where the subscript 00 denotes that the first average is to be taken with respect to the full, (λ, λ)-independent action, whereas the subscript "gh." denotes that the second average is to be taken with respect to the zero-order ghost action I 0 .Diagrammatically, for each n ≥ 1, the ghost average receives contributions from a single ghost loop, depicted in Fig. 10.
In coordinate space, suppressing the color structure, the diagram reads or, equivalently, where δ(0) is a Dirac delta in coordinate space, FIG. 10: Loop contributing to the ghost average in (A10) (example for the case n = 6).The dashed line is the (λ, λ) zero-order propagator.
Therefore, for n ≥ 1, where Ω n (ξ) is the matrix product of n factors of Ω(ξ) and the trace is taken over the color indices.
In dimensional regularization, the integral in (A13) vanishes [58].It follows that OI n 1 0,conn.= 0 for every n ≥ 1, so that, going back to (A9), where to obtain O 00 we have integrated out the free ghost action I 0 from O 0 .What (A15) means is that the per- turbative corrections to the vacuum expectation value O due to the determinant det(Λ(ξ)) vanish in dimensional regularization.Therefore, the vacuum expectation value of any operator O in the full theory can be computed by setting det(Λ(ξ)) = 1 in its dimensionally regularized partition function.
One may have noticed that our proof-aside from dimensional regularization-relies exclusively on the fact that Λ(ξ) is equal to the unit matrix to lowest order in perturbation theory.The question arises, then, whether the proof is general enough to apply to the determinant of any such matrix.The answer is that, in general, it does not.Indeed, setting δ(0) = 0 in dimensional regularization is allowed if and only if the calculations can be carried out without spoiling the symmetries of the theory.
While Lorentz invariance is clearly preserved by the action in (A3), showing that the latter does not violate the BRST invariance of the full action of the theory requires us to extend the symmetry to the ghost fields λ and λ.Indeed, a straightforward calculation starting from (11) and the definition of with so that the ghosts must have non-vanishing BRST transformations if (A3) is to be invariant.Since the BRST transformation does not act on the anti-ghost index of Λ ab (ξ), it is reasonable to define where sλ a is chosen so that s(Λλ) = 0. (A3)-and the full action of the theory together with it-is invariant with respect to these extended BRST transformations.The nilpotency of the extended BRST operator is then easily proved by observing that s 2 Λ ab (ξ) = 0-which holds thanks to the nilpotency of s on the fields ξ and c-implies that is, sΨ = −Ψ 2 .When plugged into (A18), the latter ensures that s 2 λ a = s 2 λ a = 0.

Appendix B: Unity at work
As shown in the main text, a one-loop computation shows that the σ-action with σ → σ 0 + δσ leads to where the mass term will be cancelled by the self energy contributions originating from the extra terms in the action, when using the gap equation ∂V ∂σ = 0 (cf.(43)).With the "extra" part, we mean the diagrams that are generated by the extra vertices arising from the δσ-part of the σ-Lagrangian, meaning the δσA 2 -and A 4 -vertices, see (31), (32).This will actually be a more generic feature of the dynamical model.More precisely, the tree level mass will always cancel against certain contributions coming from the extra diagrams.
Let us now try to show this by making use of the "power of unity".We will work with a simplified notation, to make things clear.But the general argument readily applies to the case under study.
Let us consider a theory, Theory 1, with a partition function [dA]e −S(A) , but we could also add a unity to get, Theory 2 with a partition function We understand that this is still the same theory as the (exact) Gaussian integral over σ yields a unity 16 .This means that (connected) correlation functions remain unchanged, A We have temporarily introduced the (loop counting) factor .Clearly, by identifying order per order in , the equivalence between A • • • A 1 and A • • • A 2 will also hold order per order in perturbation theory, seen as formal power series in .
Let us illustrate it explicitly: the extra vertices introduced by the unity, ( 31) and (32) which is indeed equal to zero.
Next, we consider the possibility that σ develops a vacuum expectation value, namely with per definition δσ = 0.In this case, the partition function will read The first exponential will still yield a unity, the second consists of a constant and will play no role in correlation functions (but it does in the quantum effective action).The term linear in δσ we dropped, as this will cancel order per order 17 With "A" a shorthand for all original fields.by removing all δσ-tadpole graphs, which is equivalent to δσ ≡ 0, or, extremization of the effective action w.r.t.σ, that is, the gap equation, see e.g.[28].Finally, there is also an effective mass term for the A-field that we included in the original action, that is, with a "dynamically massive" A-field.
At one-loop, we would now get These two diagrams are completely similar to those in the last line of (40), when resummed into the inverse propagator they will annihilate the tree level mass term upon using the gap equation ∂V ∂σ = 0, as it follows from direct computation also here.The first diagram in the first line is nonzero, but it will cancel against other (non-written) tadpole contributions.Notice that we did not explicitly use the underlying unity at this point.
We can reconsider the 3 diagrams in the first line of (B7) also from the unity viewpoint though.In that case, the 3 diagrams cancel against each other, just as before.However, this means we have removed one, but not all, tadpole contributions, meaning that there will be a net tadpole correction to AA 1PI extra , which actually corresponds to minus the first diagram, upon amputating the external A-legs of course.But due to the gap equation, this is nothing else than the tree level mass, up to the sign.
It is clear this observation will continue to hold through at any order n: we can always keep the necessary tadpole (sub)graphs to get all diagrams that make up the unity order per order (the sum of which diagrams will then lead to a zero), the remaining contribution will be the (usually never written) "counterterm" that eliminates the tadpoles, up to the sign.As at one loop, this will exactly kill the tree level mass due to the gap equation.Overall, we will thus be left by the same diagrams of the Curci-Ferrari model, up to the tree level mass.Notice though that the mass running in the loops will still need an appropriate re-expansion up to the considered order, as the actual mass is defined from m 2 ≡ −µ ε σ/ζ(g 2 ), see (25).This will lead to further differences with the Curci-Ferrari case.In (81), we constructed an RG improved potential up to leading log (LL) order.The solution (86) corresponds to a genuine stationary point as ∂V LL ∂m = 0. Upon closer inspection we notice it actually corresponds to a minimum of Re(V ), and a maximum of Im(V ).We kept however using (86) as a solution to the gap equation in the remainder of the paper 18 .Before discussing the validity of this approximation, we note that we chose to minimize V LL (m(µ)) w.r.t. to the variable m(µ).Up to the considered order, as noted below (85), one can reinterpret it in terms of V (m(µ = m)) and 18 There is actually also the solution at m 2 = 0.However, this is not to be trusted.To keep the expansion under control, we should then assume µ ∼ 0, implying g 2 ∼ ∞ in the MS scheme, thence invalidating the expansion.The relevant expansion parameter at µ 0 = 1 GeV, λ MS (µ 0 ) = 0.316, is sufficiently small to trust the used expansion more or less, with moreover m sol./µ 0 close to 1.
thence minimize w.r.t.m(m) rather than m(µ).We did not do so, as our eventual interest lies in the RG flow of the propagators in the DIS scheme, with initial conditions set at µ 0 .This would require, after having minimized w.r.t.m(m), to flow this MS solution from µ = m to µ = µ 0 using the MS anomalous dimensions, followed by converting it to the desired initial value in the DIS scheme via the conversion (87).As we do not expect the MS scheme to be trustworthy at lower scales, we intentionally first set µ = µ 0 and then directly extract an estimate for the MS initial conditions at the chosen, high enough, scale µ 0 , followed by the conversion to DIS.
In passing, we notice here that for the LL minimum solution (86), µ d dµ m sol.,MS = 0 up to the considered order.Although it might look strange to find a RG invariant solution of the LL gap equation, this will be always the case in the LL case, which after all is based on the zeroth order ("classical") potential, dressed with logs.This common feature can be easily checked using the related formalisms of [40,41] as for single scale theory, the LL potential will always have the structure of (85) with a solution similar to (86).This no-running of the solution is thus an artefact of the LL approximation.
This being said, we can now have a look at what happens when we include the next-to-leading-log (NLL) corrections.As we know the first term of the NLL series, we can use this to our advantage.Indeed, generalizing (81), we may write where w 0 = − 113N 288π 2 , as it follows from (29).The function F(u) was already determined in (84), and upon imposing (80) to next-to-leading order, we get which solves uniquely to As a check on this result, re-expanding (C1) up to order g 2 gives the correct terms upon comparison with the twoloop MS effective potential as it was computed first in [18], later on verified in [59], up to the constant term of the twoloop piece of course, which is the first of the next-to-nextto-leading log order terms.Focussing on the SU(3) case and setting as before λ MS (µ 0 ) = 0.316 at µ = µ 0 = 1 GeV, we come to Fig. 11, where both the LL and NLL potential are shown.We see that going to next-to-leading log order in the RG improved expansion does shift the location of the minimum to around 0.8 GeV, but it does land in the region where the improved NLL potential is real-valued.For the record, we refrain from using this minimum in Section III.V, as we only solved the two-point functions RG flow at leading order.Actually, the potential is rather flat in the region of interest, so let us see what happens in terms of a variable µ, which will also allow to verify if the expected UV RG asymptotics is reached.We first convert the fit value λ MS (µ 0 ) = 0.316 into the corresponding Λ MS -value through inversion of the two-loop ex-pression g 2 (µ) = 1 12 to −∞ and eventually into the complex region for too small m, making the quantity β 0 g 2 (µ) ln m 2 µ 2 too large.In any case, log RG resummations capture information from all orders, whilst the gap equation imposes a further constraint between terms of different orders.It would be interesting to investigate more sophisticated RG improvements of the effective potential, so that for example the RG flow of the solution is strictly consistent with the anomalous dimension of the mass at the chosen order, for any value of the RG scale, but this falls beyond the scope of the current paper.This might also further complicate the numerics (stability) of the fitting procedure.We plan to come back to this in future work.

( 2 )
Y M are tacitly assumed but not spelled out.

λFIG. 1 :
FIG. 1: Solutions of the gap equation in the DIS and MSschemes at the renormalization scale µ 0 = 1 GeV.
I: Parameters obtained by fitting the lattice data to the DIS scheme RG-improved gluon and ghost form factors, for different initial renormalization scales µ 0 and renormalization parameters ∆Z.The cells contain the values (λ 0 , m 0 [GeV]) of the fitted initial DIS coupling constant λ 0 = λ(µ 0 ) and of the initial DIS gluon mass parameter m 0 = m 0 (µ 0 ) computed by solving the gap equation.
µ = µ 0 , yielding Λ MS = 0.630 GeV.Feeding this and (C4) back into V NLL (m) allows one to show the potential for various values of µ, as shown in Fig.12.We observe a clear, real-valued minimum which lowers for growing µ.At last, solving for the minimum leads to Fig.13, shown together with a fitted (over the interval µ = [5, 10] GeV) consistent with the expected UV RG behaviour in terms of ∂ ln m/∂ ln µ = (γ 0 /2)g 2 +. ... For too low values of µ, i.e. too close to the MS Landau pole at µ = Λ MS , we should not trust the results anymore.The Landau pole is eventually also what pushes the potentials in FIG.