Yang-Mills propagators in linear covariant gauges from Nielsen identities

We calculate gluon and ghost propagators in Yang-Mills theory in linear covariant gauges. To that end, we utilize Nielsen identities with Landau gauge propagators and vertices as the starting point. We present and discuss numerical results for the gluon and ghost propagators for values of the gauge parameter $0<\xi \le 5$. Extrapolating the propagators to $\xi \to \infty $ we find the expected qualitative behavior. We provide arguments that our results are quantitatively reliable at least for values $\xi\lesssim 1/2$ of the gauge fixing parameter. It is shown that the correlation functions, and in particular the ghost propagator, change significantly with increasing gauge parameter. In turn, the ghost-gluon running coupling as well as the position of the zero crossing of the Schwinger function of the gluon propagator remain within the uncertainties of our calculation unchanged.


I. INTRODUCTION
In the past decades, functional approaches such as Dyson-Schwinger equations (DSEs) or functional renormalization group (fRG) equations have very successfully contributed to understanding many phenomena in quantum chromodynamics (QCD), ranging from the hadron resonance spectrum to the phase structure of QCD at nonvanishing temperatures and densities. The majority of the respective investigations have been carried out in the Landau gauge due to the technical as well as conceptual advantages this gauge provides. Evidently, obtaining via such approaches results for physical observables requires truncations to the full hierarchy of coupled functional equations, typically chosen to be of a given order in a systematic approximation scheme such as the vertex expansion. This calls for checks of the systematic errors of the respective results. For example, the gauge independence of the computed observables would be a very powerful self-consistency check. Although demonstrating generic gauge independence is likely beyond reach within functional approaches, the test of a reasonably accurate independence of gauge-invariant quantities when varying the gauge parameter within a given class of gauges would provide a convincing (but also costly) verification of the employed truncations. In a first step towards such a self-consistency check, we study in this work the propagators of elementary Yang-Mills fields in the linear covariant gauges thus extending [1]. In addition, our investigation may serve for corroborating the current state of the art of functional studies * martin.napetschnig@edu.uni-graz.at † reinhard.alkofer@uni-graz.at ‡ markus.huber@physik.jlug.de § j.pawlowski@thphys.uni-heidelberg.de in the Landau gauge, for the respective recent Landau gauge DSE results for propagators and vertex functions see, e.g., [2,3], and for recent quantitative fRG results [4,5]. These Yang-Mills results within the Landau gauge, in particular for the ghost and gluon propagators, match quantitatively the respective available lattice results, see, e.g., [6][7][8][9]. Detailed discussions of Yang-Mills correlation functions in the Landau gauge as well as more results on them can be found, for instance, in Refs. [2][3][4][5][10][11][12][13][14][15][16][17][18][19][20] and references therein. Herein, we extend previous studies [1,21] and compute ghost and gluon propagators from Nielsen identities (NIs). Further results on linear covariant gauges from functional methods can be found, e.g., in [1,[21][22][23][24][25][26][27], from the (refined) Gribov-Zwanziger framework in [28][29][30][31], from variational methods in [32][33][34], and from lattice methods in [35][36][37][38][39]; see also the respective part in the recent review [3] and references therein. The NIs describe the dependence of correlation functions on the gauge fixing parameter in terms of a differential equation of the effective action w.r.t. the gauge fixing parameter, see, e.g., [21,22,27]. The resulting equations for the correlations functions can be integrated from the Landau gauge to any linear covariant gauge, and we are going to report on an investigation in which we used the quantitative DSE results for Landau gauge correlation functions from Ref. [2] as a starting point for such an integration of a set of coupled differential equations. Moreover, NIs have also been derived for other families of gauges. Often, these families go under the name interpolating gauges. Examples include interpolating gauges between the Landau gauge and the Coulomb gauge [40,41], the Landau gauge and the maximally Abelian gauge [25,42,43], the linear covariant gauges, the Coulomb gauge and the maximally Abelian gauge [44,45], and the linear covariant gauges, the maximally Abelian gauge and the Curci-Ferrari gauge [46]. This article is structured as follows. In the next section, we introduce the correlation functions and their functional equations. In Sec. III, the setup including truncations is presented and discussed. Sec. IV contains the results, and we conclude in Sec. V with a summary. Several appendices contain technicalities including discussions of the RG and UV properties of the equations and the model parameter dependence of the solution.

II. CORRELATION FUNCTIONS AND THEIR NIELSEN IDENTITIES
As usual in functional approaches it is assumed that a Wick rotation to Euclidean space has been performed. The Lagrangian density of Yang-Mills theory in linear covariant gauges is then given by with The fields are the gluon field A a µ , the ghost field c a , the anti-ghost field c a and the Nakanishi-Lautrup field b a . The Nakanishi-Lautrup fields are introduced via the Becchi-Rouet-Stora-Tyutin (BRST) transformation [47,48], denoted by s, and given by where is the covariant derivative in the adjoint representation. Then, the gauge fixing part of the Lagrangian reads explicitly where χ = s ξ was introduced as the BRST transformation of the gauge fixing parameter ξ. The BRST transformations are nonlinear for the gluon and ghost fields. We are going to work with the Batalin-Vilkovsky (BV) or antifield formalism, see, e.g., [49][50][51][52], and introduce sources, called antifields, for them, In L BV , additional vertices appear that contain an antifield A * or c * and are hence called antifield vertices. In the Landau gauge, it is sufficient to consider the completely transverse part of the dressing functions as they form a closed system [5,12], and we split the gluon propagator into a transverse and longitudinal part written as where P µν (p) = g µν − p µ p ν /p 2 is the transverse projection operator. The Slavnov-Taylor identity (STI) for the gluon propagator enforces its longitudinal dressing function to be constant, Z L (p 2 ) = ξ, i.e., all quantum corrections are transverse. Correspondingly, in the Landau gauge, ξ = 0, the gluon propagator is proportional to a transverse projection operator. For ξ > 0, this is no longer true and also the trivial longitudinal part enters.
The ghost has only one dressing function, The ghost-gluon vertex can be conveniently split into a transverse and a longitudinal part, We use a compact notation where the subscripts denote the fields with indices corresponding to the momentum arguments. The tree-level expression is, respectively, D Acc,T (k; p, q) = 1 and D Acc,L (k; p, q) = p·k/k 2 . For the three-gluon vertex we only use a dressed tree-level tensor, thereby neglecting components which are sub-leading [53], The NI encodes the dependence of the effective action Γ on the gauge fixing parameter ξ and reads The right hand side of (11) follows from the BRST invariance of the effective action. For the final form we differentiate (11) with respect to χ and set χ = 0. This leads to a master equation for the NIs: NIs for correlation functions can now be obtained by applying further field derivatives to (12). For related work on Nielsen identities see [54,55]. For the propagators, the NIs in the form (12) used here were derived in Ref. [22], where also a perturbative oneloop analysis was done, and we refer to this reference for further details. The NIs for the ghost and gluon propagators read in Euclidean momentum space: where we have suppressed the color indices. For the numerical solution of (13) we use approximations for the vertices Γ cχA , Γ cχc * and Γ AχA * which will be discussed in detail in Sec. III.
For the ghost-gluon, the three-gluon and the four-gluon vertices the NIs read

III. TRUNCATION AND INPUT
The NIs are exact functional equations, and we cannot solve them without approximations. For the propagator equations we need Γ cχAµ (p, 0, −p), Γ cχc * (p, 0, −p) and Γ AµχA * ν (p, 0, −p). We follow Ref. [21] and calculate them from the first order in a skeleton expansion, also called dressed-loop expansion, shown in Fig. 1. An additional approximation in Ref. [21] was the use of bare vertices. In this work, however, we keep the vertices dressed. The resulting expressions K (i) arising on the right-hand side of Eq. (13) correspond to the loop diagrams i = a, b, c, d, e in Fig. 1 and are provided in Appendix A. Using these expressions and switching from the two-point functions Γ AA and Γc c to the dressing function Z and G, respectively, the propagator NIs can be written as where the i = e contribution has been split into two terms for convenience, cf., Eqs. (A1). From the loop diagrams, we can directly confirm that the NIs transform correctly under a change of the renormalization scale µ 1 → µ 2 . The NIs only depend on renormalized correlation functions and couplings O i (p, µ) that schematically change as where Z Oi (Λ, µ) is the renormalization factor of the respective correlation function or coupling O i . The explicit forms are deferred to Appendix B.
We exemplify the consistency of the left-and right-hand sides of the NIs under a change of the renormalization scale with diagram (a). Under a change of the renormalization scale µ → ν, diagram (a) transforms, upon using B2, as This is exactly as the ξ derivative on the left-hand side of Eq. (15) transforms. For all other diagrams this analysis can be repeated and leads to the same result. As starting value for integrating the NIs we use the Landau gauge, ξ = 0, for which we have results for the propagators and vertices. However, we also the need the vertices for ξ > 0. Based on the fact that the ξ dependence for the propagators found on the lattice [35][36][37][38][39] is small, we adopt as working assumption that the ghost-gluon and three-gluon vertices deviate only little as well and use results for ξ = 0 for all ξ. Only in the UV we accommodate the correct ξ dependence by modifying the anomalous running accordingly. In addition, we approximate their longitudinal parts with the transversely projected ones due to the lack of concrete results for the former. A final approximation consists in taking only a single kinematic configuration for each vertex. For the three-gluon vertex, this is a good approximation due to its small angular dependence [2,53,63,64]. The ghost-gluon vertex shows more angular dependence, which we neglect here, though. This is justified by the overall modest variation of the ghost-gluon vertex with respect to momenta. The Landau gauge results [2,65] used as initial values for solving the NIs are shown in Figs. 2 and 3 in comparison to lattice results. Specifically, we use a self-contained solution that possess several advantageous properties. Among them are manifest gauge covariance expressed by the good agreement of different couplings in the perturbative regime and a unique treatment of quadratic divergences, for details we refer to [2]. We also fix the overall scale from these results which which was obtained by matching the maximum of the gluon dressing function to lattice results. In the Landau gauge it is well studied that a family of different solutions can be obtained from functional equations which differ in their IR behavior [12,[66][67][68][69][70][71]. These solutions are different only in the region below 2 GeV. Most notably, the maximum in the gluon propagator changes. To explore the existence of such solutions also beyond the Landau gauge we choose two different sets of solutions. One of the solutions is the one that agrees best with lattice results and has only a shallow, hardly visible maximum in the gluon propagator. As a second choice, we take a solution with a more pronounced maximum, see Fig. 2. It remains to specify the models used for the vertices Γ A * Ac and Γ c * cc . When using bare vertices, we found that in the infrared (IR) individual loop diagrams can qualitatively modify the IR solution for the gluon propagator. This is either resolved by cancellations between individual diagrams or by the IR behavior of the antifield vertices. We explore the second option, as it is currently not clear if the first one can be realized in the truncation we use.
As Ansatz for the vertices we multiply their tree-level tensors with products of the ghost and gluon dressing functions with appropriate powers, with p 2 = (p 2 + q 2 + k 2 )/2. The denominator ensures that the vertex models are unity at p 2 = s. The antifield vertices run logarithmically like the ghost-gluon vertex, as can be checked by a perturbative one-loop analysis.
The exponents are determined such that they respect this UV behavior. As a second condition, the integrals in the NIs should be IR finite. We make the simple Ansatz where α i and β i are ξ-independent parameters. Enforcing the conditions above, we obtain β 0 > 0 is a free parameter for which we choose for convenience β 0 = 1. A test of the sensitivity of our results on this choice as well as a discussion of the function H(x) is provided in Appendix D. It can be seen that the parameter β 0 , if chosen within a reasonable range, only influences the IR and this in a quantitatively mild way. Both the ghost and the gluon NI have the form K denotes the integrals from the skeleton expansions of the vertices. The formal solution to this equation is where ξ 0 denotes the gauge fixing parameter for a known solution. The quantity K is obtained by numerically calculating the integrals which are standard one-loop twopoint integrals. As a computational framework we use CrasyDSE [72]. The integrals are logarithmically divergent and need to be renormalized. We do so by momentum subtraction and require that the dressing functions stay the same at the highest calculated momentum point. It should be noted that no quadratic divergences [73] are present. This situation has to be contrasted with the DSE or flow equation for the gluon propagator for which spurious quadratic divergences can arise, cf. Refs. [2,4,74].

IV. RESULTS
For the full nonperturbative solution, we solve the differential equations (13) up to ξ = 5 starting from the Landau gauge. The check of the self-consistency of the UV limit is deferred to Appendix C.
The results for the gluon dressing functions and the gluon propagators are shown in Fig. 4. With increasing ξ, the gluon propagator decreases. This can be seen in both the maximum of the dressing function and the IR behavior of the propagator. The difference between Landau gauge and ξ = 0.5 is not drastic and compatible with lattice results [37]. We find at 0 and 1 GeV that the gluon propagator goes down by 8% and 7%, respectively, for ξ = 0.5. For the lattice results these ratios are given in Ref. [37] as approximately 10% and 5% with errors of a few percentage points each. It should be noted, though, that the IR behavior of the gluon propagator depends on the employed models for the antifield vertices. All individual diagrams in the gluon NI are IR finite as shown in Fig.  12 of Appendix A. This comes from the IR behavior of the anti-field vertices. If we used bare vertices, IR divergences would arise that would qualitatively change the IR behavior of the gluon propagator.
The ghost propagator shows for ξ > 0 the already known logarithmic IR suppression [1,21], see Fig. 5. This behavior results from diagram (a), see Fig. 12. For small ξ, the deviation from the Landau gauge is not very large. This agrees with lattice results which do not see a change  in the ghost propagator up to ξ = 0.3 and above approximately 500 MeV, which was the lowest accessible momentum value [38,39]. In the continuum results displayed here, evaluated down to 0.01 MeV, we do, however, see deviations from the Landau gauge behavior below 500 MeV also for small values of the gauge parameter ξ. Sizeable deviations appear for higher values of the gauge fixing parameter for which also the effect on the UV behavior becomes visible.
Given the comparatively simple approximation employed for the NIs, the agreement with lattice results is very good. On the other hand, the method is very stable, and we calculated up to ξ = 5 without encountering any problems. Indeed, we can easily check that for ξ = 3 and ξ = 13/3 the correct one-loop anomalous dimensions are produced, as for these values the ghost and the gluon anomalous dimensions vanish, respectively. This is illustrated in Fig. 6 where the corresponding propagators are compared to the Landau gauge ones in the momentum region from 1 to 10 5 GeV 2 . The ghost-gluon coupling is defined via the relation where the ghost-gluon vertex dressing D Acc is evaluated at the symmetric point, and α(µ 2 ) = g 2 /(4π). One-loop universality entails that any dependence of the coupling on the gauge fixing parameter is suppressed at high momenta. Beyond one-loop, a dependence on ξ can appear, see, e.g., Refs. [75][76][77]. However, we can still assess the effect of the truncation by comparing the couplings in the perturbative regime above a few GeV. For the coupling the correct running of all quantities is important. We thus use the one-loop resummed expression for the ghost-gluon vertex, where ω = 11 N c α(s)/(12π)G 2 (s)Z(s)[D Acc (s)] 2 . We show the couplings for various values of ξ in Fig. 7.  Hereby, the scale s is chosen as 10 5 GeV 2 , phrased otherwise, the couplings are fixed at this value. As can be seen, down to 10 GeV they agree even for higher values of the gauge fixing parameters. Up to ξ = 0.5, the agreement is good down to approximately 2 GeV. To appreciate this agreement, we also show the coupling without the ghostgluon vertex dressing in Fig. 7. The agreement is worse then and the order of magnitudes is different with ξ = 0 having the largest coupling. We additionally show the one-loop vertex dressing in the inset to highlight that the vertex dressing becomes sizable already for low ξ. It would be interesting to test the ξ dependence of related quantites like the effective charge defined in [78]. Another interesting quantity is the Schwinger function ∆(t) of the gluon propagator, defined as the Fourier transformation of the momentum space propagator for vanishing spatial momentum. If the propagator violates positivity, this is reflected in the Schwinger function. Fig. 8 shows the Schwinger function for various values of the gauge fixing parameter. Up to approximately ξ = 1 the Schwinger function barely changes. In particular, the position of the zero crossing does not move. Only for higher values of ξ it moves slightly. This stability is in marked contrast to the situation for the family of solutions in the Landau gauge where the position of the zero crossing moves [2]. This can be understood as the existence of the zero crossing is related to the maximum of the gluon propagator. For linear covariant gauges, we find that the position of this maximum is basically constant in ξ. Different members of the family of solutions in the Landau gauge, on the other hand, exhibit different positions for the maxima [2,4] and hence the Schwinger function also changes. Finally, to explore the fate of the family of different solutions for correlation functions in the Landau gauge [12,16,71,[79][80][81], we also solved the NIs using a second Landau gauge solution. In this context it should be mentioned that these different solutions may correspond to different nonperturbative infrared gauge completions of the perturbative Landau gauge as discussed in [12,80]. If this conjecture is correct, this would correspond to a second gauge fixing direction in addition to ξ. Indeed, all observables in Yang-Mills theory and QCD, computed so far within this potential family of infrared completions of the Landau gauge agree within the respective error bars. A specifically relevant example in the present context of Yang-Mills theories is provided by the glueball masses, see [82]. In line with the conjecture discussed above, the obtained masses did not show any deviations within errors. In all plots in this section we used up to here the solution that is closest to lattice results. It is characterized by a very flat maximum of the gluon propagator, and a ghost dressing function that is relatively small in the deep infrared. The different Landau gauge solution used next as a starting point for the NIs is shown in Fig. 2 in comparison to the previously used solution. The second solution has a pronounced maximum in the gluon propagator and shows a clear increase of the ghost dressing function at low momenta. The results for a selection of values of ξ are shown in Fig. 9. The typical features of the Landau gauge solution type is inherited by the ξ > 0 ones. In particular, the gluon propagator has a shallow/pronounced maximum from which it follows immediately that it violates positivity. Note that such a property also leads to a spectral dimension of one in the deep IR [83]. Correspondingly, if the maximum vanished, this would imply a qualitative change of the type of solution, and it is reassuring that we do not observe that. Since a nonzero gauge fixing parameter washes out the gauge fixing condition of the Landau gauge and thus the differences between the two solutions, it is interesting to check if the two solutions approach each other for high values of ξ. To assess that, we plot the ghost dressing function and the gluon propagator at fixed momenta as a function of ξ in Fig. 10. For the ghost dressing function we see that the two solutions come closer to each other for higher values. For the gluon propagator, on the other hand, this effect is not observed. In both cases it seems plausible that for ξ → ∞ the functions vanish. This limit corresponds to removing the gauge fixing which in fact necessarily will eventually lead to a vanishing gluon propagator. Phrased otherwise, we see the expected behavior based on the general properties of the linear covariant gauges. This provides some confidence that the overall qualitative behavior of the propagators is correct for all allowable values of the gauge parameter 0 < ξ < ∞.
The distinct relative behavior of the two solutions is clearly visible in Fig. 11 which shows the ratio between the propagators at a fixed momentum point in the IR. The ratio stays constant for the gluon propagator but depends strongly on the gauge fixing parameter for the ghost propagator. Since the plotted ratio is G 1 /G 2 , this means that the ghost propagator with higher values in the IR for the Landau gauge decreases faster than the one with lower values.

V. SUMMARY
We have calculated the ghost and gluon propagators of Yang-Mills theory, respectively, quenched QCD, in the linear covariant gauges for values of the gauge fixing parameter 0 < ξ ≤ 5. The starting point has been results in the Landau gauge, ξ = 0, which were obtained in a self-contained DSE calculation. As external input we employed the nonperturbative parts of the ghost-gluon and three-gluon vertices of Landau gauge for all values of ξ. In addition, we used ansätze for the antifield vertices which contain one free parameter. We found that the solutions are not very sensitive to variations of this parameter.
In the IR, we recover the logarithmic suppression of the ghost dressing function, predicted by earlier investigations, and find an IR finite gluon propagator. The latter result, however, happens by construction based on the antifield vertex model. All results agree well, even quantitatively, with available lattice results. Compared to other methods, our setup is quite stable, even at values of ξ beyond the Feynman gauge. In particular, we recover the correct UV behavior for the propagators most convincingly seen by vanishing anomalous dimensions for the ghost and gluon dressing functions at ξ = 3 (Yennie gauge) and ξ = 13/3, respectively. We did not encounter any signs of instability up to the highest calculated value, ξ = 5. While the changes of propagators and vertices are sizable, observables are ξ-independent. This calls for respective studies of e.g. glueball masses as done in [82] for the Landau gauge. Such a study was beyond the scope of the present work. Instead, as a first step in this direction, we have discussed the ξ dependence of the ghost-gluon coupling, Fig. 7, and the zero crossing of the Schwinger function, Fig. 8. We have shown that the ξ dependence of both, the coupling as well as the Schwinger function zero crossing, are very small up to ξ = 0.5, which is highly nontrivial. Beyond ξ = 0.5, the reliability of the current approximation is successively getting worse because we do not consider the back-coupling of the ξ dependence in the vertices. Nevertheless, the observed deviations are still quite small. We also have explored the potential family of nonperturbative infrared completions, as discussed in the Landau gauge. We have tested two different starting points and obtained two corresponding sets of solutions for ξ > 0. The qualitative features, in particular violation of positivity, remain intact at least up to ξ = 5. In the limit of infinite ξ, both propagators are in agreement with the expectation that they vanish in this limit.
In the present work we have restricted ourselves to pure Yang-Mills theory. However, the inclusion of dynamical quarks is straightforward as there are no direct quark contributions in the gluon and ghost Nielsen identities [22]. Moreover, the Nielsen identity for the quark propagator has a similar structure as those for the other propagators and could be solved within a skeleton expansion. It would be also interesting to extend the current study to other covariant gauges like the maximally Abelian gauge. There, direct calculations are complicated due to its IR dominant two-loop diagrams [25,84].

ACKNOWLEDGMENTS
We thank Joannis Papavassiliou for discussions. Support by the FWF (Austrian science fund) under The loops displayed in Fig. 1 lead to the following expressions: G(q)G(p + q) p 2 q 4 (p + q) 2 (p 2 q 2 − (p · q) 2 ))D Acc T (−p; p + q, q)Γ AA * c (−q, −p, p + q), Z(p + q)G(q) p 2 q 4 (p + q) 4 (q 2 + 2 p · q)(3 p 4 + (p · q) 2 + 2 p 2 (q 2 + 3 p · q)) × C AAA (−p, p + q, −q)Γ AA * c (p + q, −p, −q), (A1e) G(q) q 4 (p + q) 4 (p 2 q 2 − (p · q) 2 )C AAA (−p, p + q, −q)Γ AA * c (p + q, −p, −q). (A1f) The last integral was split into two parts to disentangle the contribution from the transverse and longitudinal parts of the gluon propagator. We used p µ q ν r ρ Γ AAA µνρ (p, q, r) = 0 in several places to simplify the expressions. E.g., due to this only the transversely projected part of the ghost-gluon vertex appears in the second diagram. The results for the individual K (i) are shown in Fig. 12 for two different values of ξ. The model we employ for the antifield vertices depends on one parameter β 0 . For the results shown in the main part we used β 0 = 1. We tested the influence of this parameter by calculating the propagators also with β 0 = 0.5. The model function H(p 2 ) for these two choices is shown in Fig. 13. In the quantitatively relevant regime around 1 GeV the two parameter values lead only to small differences for all ξ. The rise in the UV for larger ξ comes directly from the anomalous dimension of the antifield vertices.
The propagators obtained from β 0 = 0.5 are compared to the ones from β 0 = 1 in Fig. 14. We can see that changing β 0 affects basically only the IR. Only for ξ > 4 small effects in the ghost dressing function are seen also in the midmomentum regime. We thus conclude that within the present approximation scheme the dependence on the model for the antifield vertices is of minor importance and only seen quantitatively for low momenta.