# Yang-Mills propagators in linear covariant gauges from Nielsen identities

Martin Napetschnig<EMAIL_ADDRESS>Institute of Physics,
University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz, Austria
Reinhard Alkofer<EMAIL_ADDRESS>Institute of Physics,
University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz, Austria  Markus
Q. Huber<EMAIL_ADDRESS>Institut für Theoretische Physik,
Justus-Liebig-Universität Giessen, 35392 Giessen, Germany  Jan M. Pawlowski
<EMAIL_ADDRESS>Institute for Theoretical Physics,
Universität Heidelberg, Philosophenweg 12, D-69120 Germany  ExtreMe Matter
Institute EMMI, GSI, Planckstr. 1, D-64291 Darmstadt, Germany

###### Abstract

We calculate gluon and ghost propagators in Yang-Mills theory in linear
covariant gauges. To that end, we utilise 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\leq 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.

correlation functions, Yang-Mills theory, linear covariant gauges, Nielsen
identities

###### pacs:

12.38.Aw, 14.70.Dj, 12.38.Lg

##  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 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. [10, 11, 12, 13, 14, 15, 16, 17, 18, 4, 19,
3, 2, 20, 5] and references therein.

Herein, we extend previous studies [21, 1] 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 [22, 23, 24, 25, 21, 1,
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., [22, 21, 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 [42, 25, 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

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{\text{YM}}+\mathcal{L}_{\text{gf}}\,,$ (1)

with

$\displaystyle\mathcal{L}_{\text{YM}}$
$\displaystyle=\frac{1}{4}F_{\mu\nu}^{a}F^{a}_{\mu\nu}\,,$ $\displaystyle
F^{a}_{\mu\nu}$
$\displaystyle=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}-g\,f^{abc}A_{\mu}^{b}A_{\nu}^{c}\,,$
$\displaystyle\mathcal{L}_{\text{gf}}$
$\displaystyle=s\left(\overline{c}^{a}\partial_{\mu}A_{\mu}^{a}-i\frac{\xi}{2}\overline{c}^{a}b^{a}\right)\,.$
(2)

The fields are the gluon field $A_{\mu}^{a}$, the ghost field $c^{a}$, the
anti-ghost field $\overline{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

$\displaystyle s\,A_{\mu}^{a}$ $\displaystyle=-D_{\mu}^{ab}c^{b},$
$\displaystyle s\,c^{a}$ $\displaystyle=-\frac{1}{2}g\,f^{abc}c^{b}c^{c},$
$\displaystyle s\,\bar{c}^{a}$ $\displaystyle=i\,b^{a},$ $\displaystyle
s\,b^{a}$ $\displaystyle=0\,,$ (3)

where

$\displaystyle D_{\mu}^{ab}$
$\displaystyle=\delta^{ab}\partial_{\mu}+g\,f^{abc}A_{\mu}^{c}\,,$ (4)

is the covariant derivative in the adjoint representation. Then, the gauge
fixing part of the Lagrangian reads explicitly

$\displaystyle\mathcal{L}_{\text{gf}}$
$\displaystyle=i\,b^{a}\,\partial_{\mu}A^{a}_{\mu}-\overline{c}^{a}\,\partial_{\mu}(-D^{ab}_{\mu}c^{b})-\frac{i}{2}\chi\,\overline{c}^{a}\,b^{a}\,,$
(5)

where $\chi=s\,\xi$ was introduced as the BRST transformation of the gauge
fixing parameter $\xi$.

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,

$\displaystyle\mathcal{L}_{\text{BV}}$
$\displaystyle=-A_{\mu}^{*a}s\,A^{a}_{\mu}-c^{a*}s\,c^{a}$
$\displaystyle=A_{\mu}^{*a}D^{ab}_{\mu}c^{b}+\frac{g}{2}f^{abc}c^{*a}c^{b}c^{c}.$
(6)

In $\mathcal{L}_{\text{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 [12, 5], and we
split the gluon propagator into a transverse and longitudinal part written as

$\displaystyle D^{ab}_{\mu\nu}(p)$
$\displaystyle=\delta^{ab}D_{\mu\nu}(p)=\delta^{ab}(D^{T}_{\mu\nu}(p)+D^{L}_{\mu\nu}(p))\,,$
$\displaystyle D^{T}_{\mu\nu}(p)$
$\displaystyle=P_{\mu\nu}(p)\frac{Z(p^{2})}{p^{2}}\,,$ $\displaystyle
D^{L}_{\mu\nu}(p)$
$\displaystyle=\frac{p_{\mu}p_{\nu}}{p^{2}}\frac{Z^{L}(p^{2})}{p^{2}}\,,$ (7)

where $P_{\mu\nu}(p)=g_{\mu\nu}-p_{\mu}p_{\nu}/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})=\xi$, i.e., all quantum corrections are transverse.
Correspondingly, in the Landau gauge, $\xi=0$, the gluon propagator is
proportional to a transverse projection operator. For $\xi>0$, this is no
longer true and also the trivial longitudinal part enters.

The ghost has only one dressing function,

$\displaystyle D^{ab}_{G}(p)=-\delta^{ab}\frac{G(p^{2})}{p^{2}}\,.$ (8)

The ghost-gluon vertex can be conveniently split into a transverse and a
longitudinal part,

$\displaystyle\Gamma_{A^{a}_{\mu}\bar{c}^{b}c^{c}}(k;p,q)=$
$\displaystyle\,-i\,g\,f^{abc}\Bigl{(}D^{A\bar{c}c,T}(k;p,q)P_{\mu\nu}(k)p_{\nu}$
$\displaystyle\hskip 28.45274pt+D^{A\bar{c}c,L}(k;p,q)k_{\mu}\Bigr{)}\,.$ (9)

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^{A\bar{c}c,T}(k;p,q)=1$ and $D^{A\bar{c}c,L}(k;p,q)=p\cdot
k/k^{2}$. For the three-gluon vertex we only use a dressed tree-level tensor,
thereby neglecting components which are sub-leading [53],

$\displaystyle\Gamma_{A^{a}_{\mu}A^{b}_{\nu}A^{c}_{\rho}}(p,q,r)=i\,g\,f^{abc}C^{AAA}(p,q,r)$
$\displaystyle\quad\times\Big{[}(p-q)_{\rho}g_{\mu\nu}+(q-r)_{\mu}g_{\nu\rho}+(r-p)_{\nu}g_{\mu\rho}\Big{]}\,.$
(10)

The NI encodes the dependence of the effective action $\Gamma$ on the gauge
fixing parameter $\xi$ and reads

$\displaystyle\frac{\partial\Gamma}{\partial\xi}\chi=\int
dx\left(\frac{\delta\Gamma}{\delta A_{\mu}^{a}}\frac{\delta\Gamma}{\delta
A_{\mu}^{*a}}+\frac{\delta\Gamma}{\delta c^{a}}\frac{\delta\Gamma}{\delta
c^{*a}}+i\frac{\delta\Gamma}{\delta\bar{c}^{a}}b^{a}\right)\,.$ (11)

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 $\chi$ and
set $\chi=0$. This leads to a master equation for the NIs:

$\displaystyle\left.\frac{\partial\Gamma}{\partial\xi}\right|_{\chi=0}=\int
dx\left.\left(\frac{\partial\delta\Gamma}{\partial\chi\delta
A_{\mu}^{a}}\frac{\delta\Gamma}{\delta
A^{*a}_{\mu}}+\frac{\delta\Gamma}{\delta
A_{\mu}^{a}}\frac{\partial\delta\Gamma}{\partial\chi\delta
A^{*a}_{\mu}}\right.\right.$
$\displaystyle\left.\left.+\frac{\partial\delta\Gamma}{\partial\chi\delta
c^{a}}\frac{\delta\Gamma}{\delta
c^{{}^{*}a}}-\frac{\delta\Gamma}{\delta{c^{a}}}\frac{\partial\delta\Gamma}{\partial\chi\delta
c^{*a}}+ib^{a}\frac{\partial\delta\Gamma}{\partial\chi\delta\bar{c}^{a}}\right)\right|_{\chi=0}\,.$
(12)

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 one-loop 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:

$\displaystyle\partial_{\xi}\Gamma_{\bar{c}c}(p^{2})=$
$\displaystyle\,\bigg{[}i\frac{p_{\mu}}{p^{2}}\Gamma_{\overline{c}\chi
A_{\mu}}(p,0,-p)$ (13a) $\displaystyle\hskip 11.38092pt+\Gamma_{c\chi
c^{*}}(p,0,-p)\bigg{]}\Gamma_{\bar{c}c}(p^{2})\,,$
$\displaystyle\partial_{\xi}\Gamma_{A_{\mu}A_{\nu}}(p^{2})=$
$\displaystyle\,2\Gamma_{A^{*}_{\rho}\chi
A_{\mu}}(p,0,-p)\Gamma_{A_{\nu}A_{\rho}}(p^{2})\,,$ (13b)

where we have suppressed the color indices. For the numerical solution of (13)
we use approximations for the vertices $\Gamma_{\overline{c}\chi A}$,
$\Gamma_{c\chi c^{*}}$ and $\Gamma_{A\chi 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

$\displaystyle\partial_{\xi}\>\Gamma_{A_{\mu}^{a}\overline{c}^{b}c^{c}}$
$\displaystyle=-i\,\frac{p_{\rho}}{p^{2}}\Gamma_{\chi\overline{c}^{b}A^{d}_{\rho}A^{a}_{\mu}}\Gamma_{\overline{c}^{d}c^{b}}+\Gamma_{\chi\overline{c}^{b}A^{d}_{\rho}}\Gamma_{c^{c}A^{*d}_{\rho}A^{a}_{\mu}}+\Gamma_{c^{c}\chi
c^{*d}}\Gamma_{A^{a}_{\mu}\overline{c}^{b}c^{d}}$
$\displaystyle\quad+\Gamma_{\overline{c}^{b}c^{d}}\Gamma_{A^{a}_{\mu}c^{c}\chi
c^{*d}}+\Gamma_{A^{a}_{\mu}A^{d}_{\rho}}\Gamma_{\overline{c}^{b}c^{c}\chi
A^{*d}_{\rho}}+\Gamma_{A^{d}_{\rho}\overline{c}^{b}c^{c}}\Gamma_{\chi
A^{a}_{\mu}A^{*d}_{\rho}}\,,$
$\displaystyle\partial_{\xi}\>\Gamma_{A_{\mu}^{a}A_{\nu}^{b}A_{\rho}^{c}}$
$\displaystyle=\Gamma_{\chi
A_{\sigma}^{*d}A_{\mu}^{a}A_{\nu}^{b}}\Gamma_{A_{\rho}^{c}A_{\sigma}^{d}}+\Gamma_{\chi
A_{\sigma}^{*d}A_{\mu}^{a}}\Gamma_{A_{\nu}^{b}A_{\rho}^{c}A_{\sigma}^{d}}+\text{permutations}\,,$
$\displaystyle\partial_{\xi}\>\Gamma_{A_{\mu}^{a}A_{\nu}^{b}A_{\rho}^{c}A_{\sigma}^{d}}$
$\displaystyle=\Gamma_{\chi
A^{*e}_{\tau}A_{\mu}^{a}A_{\nu}^{b}A_{\rho}^{c}}\Gamma_{A_{\sigma}^{d}A_{\tau}^{e}}+\Gamma_{\chi
A^{*e}_{\tau}A_{\mu}^{a}A_{\nu}^{b}}\Gamma_{A_{\rho}^{c}A_{\sigma}^{d}A_{\tau}^{e}}+\Gamma_{\chi
A^{*e}_{\tau}A_{\mu}^{a}}\Gamma_{A_{\nu}^{b}A_{\rho}^{c}A_{\sigma}^{d}A_{\tau}^{e}}+\text{permutations}\,.$
(14)

##  III Truncation and input

(a) (b)

(c) (b)

(d) (e)

Figure 1: Skeleton expansions for $\Gamma_{\overline{c}\chi A_{\mu}}$,
$\Gamma_{c\chi c^{*}}$, and $\Gamma_{A\chi A^{*}}$ (top to bottom). Wiggly
lines denote gluons and dotted ones ghosts. Wiggly-dashed ones represent a
mixed gluon-Nakanishi Lautrup field propagator. The antifields $A^{*}$ and
$c^{*}$ as well as $\chi$ are indicated explicitly. Propagators are all
dressed, small/large dots represent bare/dressed vertices.

The NIs are exact functional equations, and we cannot solve them without
approximations. For the propagator equations we need $\Gamma_{\overline{c}\chi
A_{\mu}}(p,0,-p)$, $\Gamma_{c\chi c^{*}}(p,0,-p)$ and $\Gamma_{A_{\mu}\chi
A^{*}_{\nu}}(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
$\Gamma_{AA}$ and $\Gamma_{\bar{c}c}$ to the dressing function $Z$ and $G$,
respectively, the propagator NIs can be written as

$\displaystyle\partial_{\xi}\ln G(p^{2})$
$\displaystyle=K^{(a)}+K^{(b)}+K^{(c)},$ $\displaystyle\partial_{\xi}\ln
Z(p^{2})$ $\displaystyle=K^{(d)}+K^{(e_{1})}+K^{(e_{2})}\,,$ (15)

where the $i=e$ contribution has been split into two terms for convenience,
cf., Eqs. (25).

Figure 2: Ghost and gluon dressing functions (left) and gluon propagator
(right) from lattice [56] and DSE calculations [2] in the Landau gauge.

Figure 3: Ghost-gluon (left) and three-gluon vertices (right) in Landau gauge
from lattice [57, 58, 59] and DSE calculations [2] in the Landau gauge. For
more lattice results on the three-gluon vertex see, e.g., [60, 61, 62].

Figure 4: The gluon dressing function (left) and propagator (right) for
various $\xi$ including the starting point $\xi=0$.

From the loop diagrams, we can directly confirm that the NIs transform
correctly under a change of the renormalization scale $\mu_{1}\to\mu_{2}$. The
NIs only depend on renormalized correlation functions and couplings ${\cal
O}_{i}(p,\mu)$ that schematically change as

$\displaystyle{\cal O}_{i}(p,\mu_{1})\to\frac{Z_{{\cal
O}_{i}}(\Lambda,\mu_{1})}{Z_{{\cal O}_{i}}(\Lambda,\mu_{2})}{\cal
O}_{i}(p,\mu_{2})\,,$ (16)

where $Z_{{\cal O}_{i}}(\Lambda,\mu)$ is the renormalization factor of the
respective correlation function or coupling ${\cal 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 $\mu\rightarrow\nu$, diagram (a) transforms, upon
using 27, as

$\displaystyle\frac{Z_{g}(\Lambda,\nu)^{2}\widetilde{Z}_{3}(\Lambda,\nu)^{2}\widetilde{Z}_{1}(\Lambda,\mu)^{2}}{Z_{g}(\Lambda,\mu)^{2}\widetilde{Z}_{3}(\Lambda,\mu)^{2}\widetilde{Z}_{1}(\Lambda,\nu)^{2}}=\frac{Z_{3}(\Lambda,\mu)}{Z_{3}(\Lambda,\nu)}.$
(17)

This is exactly as the $\xi$ 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, $\xi=0$,
for which we have results for the propagators and vertices. However, we also
the need the vertices for $\xi>0$. Based on the fact that the $\xi$ 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 $\xi=0$ for all $\xi$. Only in
the UV we accommodate the correct $\xi$ 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 [63, 53, 2, 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
[66, 67, 68, 69, 70, 12, 71]. These solutions are different only in the region
below $2\,\text{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 $\Gamma_{A^{*}Ac}$ and
$\Gamma_{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.

Figure 5: The ghost dressing function for various $\xi$ including the starting
point $\xi=0$. Figure 6: The UV behavior of the ghost and gluon dressing
functions for $\xi=3$ and $\xi=13/3$, respectively.

As Ansatz for the vertices we multiply their tree-level tensors with products
of the ghost and gluon dressing functions with appropriate powers,

$\displaystyle\Gamma_{A^{*a}_{\mu}A_{\nu}^{b}c^{c}}(p,q,k)$
$\displaystyle=g\,f^{abc}\,H(\overline{p}^{2})\,g_{\mu\nu}\,,$
$\displaystyle\Gamma_{c^{*a}c^{b}c^{c}}(p,q,k)$
$\displaystyle=-g\,f^{abc}\,H(\overline{p}^{2})\,,$ $\displaystyle H(x)$
$\displaystyle=\frac{G(x)^{\alpha}Z(x)^{\beta}}{G(s)^{\alpha}Z(s)^{\beta}}\,,$
(18)

with $\overline{p}^{2}=(p^{2}+q^{2}+k^{2})/2$. The denominator ensures that
the vertex models are unity at $\overline{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

$\displaystyle\alpha=\alpha_{0}+\alpha_{1}\xi\,,\qquad\beta=\beta_{0}+\beta_{1}\xi\,,$
(19)

where $\alpha_{i}$ and $\beta_{i}$ are $\xi$-independent parameters. Enforcing
the conditions above, we obtain

$\displaystyle\alpha_{0}=-\frac{26\beta_{0}}{9}\,,\quad\alpha_{1}=-\frac{9-4\beta_{0}}{6}\,,\quad\beta_{1}=\frac{9-4\beta_{0}}{12}\,.$
(20)

$\beta_{0}>0$ is a free parameter for which we choose for convenience
$\beta_{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 $\beta_{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

$\displaystyle\frac{\partial\ln M}{\partial\xi}=K\,.$ (21)

$K$ denotes the integrals from the skeleton expansions of the vertices. The
formal solution to this equation is

$\displaystyle M(\xi)=M(\xi_{0})e^{\int_{\xi_{0}}^{\xi}d\xi K}\,,$ (22)

where $\xi_{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 two-point 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. [4, 2, 74].

##  IV Results

Figure 7: Left: The coupling for various values of $\xi$. Right: The coupling
without the vertex dressing. The inset shows the one-loop expressions for the
ghost-gluon vertex, Eq. (24), used for the couplings.

For the full nonperturbative solution, we solve the differential equations
(13) up to $\xi=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 $\xi$, 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 $\xi=0.5$ is not
drastic and compatible with lattice results [37]. We find at $0$ and
$1\,\text{GeV}$ that the gluon propagator goes down by $8\%$ and $7\%$,
respectively, for $\xi=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 $\xi>0$ the already known logarithmic IR
suppression [21, 1], see Fig. 5. This behavior results from diagram (a), see
Fig. 12. For small $\xi$, 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 $\xi=0.3$ and above approximately $500\,\text{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 $\xi$. 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 $\xi=5$ without encountering any
problems. Indeed, we can easily check that for $\xi=3$ and $\xi=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}$ GeV2.

The ghost-gluon coupling is defined via the relation

$\displaystyle\alpha(p^{2})=\alpha(\mu^{2})G^{2}(p^{2})Z(p^{2})[D^{A\bar{c}c}(p^{2})]^{2}\,,$
(23)

where the ghost-gluon vertex dressing $D^{A\bar{c}c}$ is evaluated at the
symmetric point, and $\alpha(\mu^{2})=g^{2}/(4\pi)$. 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 $\xi$ 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,

$\displaystyle
D^{A\bar{c}c,1l}(p^{2})=D^{A\bar{c}c}(s)\left(1+\omega\ln\frac{p^{2}}{s}\right)^{-3\xi/22}\,,$
(24)

where $\omega=11\,N_{c}\,\alpha(s)/(12\pi)G^{2}(s)Z(s)[D^{A\bar{c}c}(s)]^{2}$.
We show the couplings for various values of $\xi$ in Fig. 7. Hereby, the scale
$s$ is chosen as $10^{5}\,\text{GeV}^{2}$, phrased otherwise, the couplings
are fixed at this value. As can be seen, down to $\lesssim$ 10 GeV they agree
even for higher values of the gauge fixing parameters. Up to $\xi=0.5$, the
agreement is good down to approximately 2 GeV. To appreciate this agreement,
we also show the coupling without the ghost-gluon vertex dressing in Fig. 7.
The agreement is worse then and the order of magnitudes is different with
$\xi=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 $\xi$. It would be interesting to test the $\xi$ dependence of
related quantites like the effective charge defined in [78].

Figure 8: The absolute value of the Schwinger function for various values of
the gauge fixing parameter. The Schwinger function is negative for $t\gtrsim
1.3\,\text{fm}$.

Figure 9: Ghost dressing function (left) and gluon propagator (right) for two
different decoupling solutions as starting points at various value of $\xi$.
Figure 10: The gauge parameter dependence of the ghost dressing function
(dashed line) and gluon propagator (solid line) at fixed momenta. The lowest
calculated value is for $1/\xi=1/5$. Figure 11: The ratio of the propagators
for two solutions. For the ghost propagator the lowest calculated momentum is
used, for the gluon propagator zero momentum.

Another interesting quantity is the Schwinger function $\Delta(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 $\xi=1$ the Schwinger function barely changes. In particular,
the position of the zero crossing does not move. Only for higher values of
$\xi$ 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 $\xi$. Different members of the family of solutions in the Landau
gauge, on the other hand, exhibit different positions for the maxima [4, 2]
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 [79, 12, 71, 80, 16, 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 $\xi$. 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 $\xi$ are shown in Fig. 9. The typical features
of the Landau gauge solution type is inherited by the $\xi>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 $\xi$. To assess that, we plot the ghost dressing function and the gluon
propagator at fixed momenta as a function of $\xi$ 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 $\xi\rightarrow\infty$ 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<\xi<\infty$.

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<\xi\leq 5$. The starting point has been results in
the Landau gauge, $\xi=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 $\xi$.
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 $\xi$
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 $\xi=3$ (Yennie gauge) and
$\xi=13/3$, respectively. We did not encounter any signs of instability up to
the highest calculated value, $\xi=5$.

While the changes of propagators and vertices are sizable, observables are
$\xi$-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 $\xi$ dependence of the ghost-gluon coupling, Fig. 7, and the zero
crossing of the Schwinger function, Fig. 8. We have shown that the $\xi$
dependence of both, the coupling as well as the Schwinger function zero
crossing, are very small up to $\xi=0.5$, which is highly nontrivial. Beyond
$\xi=0.5$, the reliability of the current approximation is successively
getting worse because we do not consider the back-coupling of the $\xi$
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 $\xi>0$.
The qualitative features, in particular violation of positivity, remain intact
at least up to $\xi=5$. In the limit of infinite $\xi$, 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 Contract No. P27380-N27 is gratefully acknowledged. This
work is supported by EMMI and the BMBF grant 05P18VHFCA. It is part of and
supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT) and the
DFG under Germany’s Excellence Strategy EXC - 2181/1 - 390900948 (the
Heidelberg Excellence Cluster STRUCTURES).

##  Appendix A Diagrams in the Nielsen identities of the propagators

The loops displayed in Fig. 1 lead to the following expressions:

$\displaystyle K^{(a)}$
$\displaystyle=-\frac{N_{c}g^{2}}{2}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{G(q)G(p+q)}{q^{2}(p+q)^{2}}D_{L}^{A\overline{c}c}(-p;p+q,-q)D_{L}^{A\overline{c}c}(q;p,-p-q),$
(25a) $\displaystyle K^{(b)}$
$\displaystyle=-\frac{N_{c}g^{2}}{2}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{G(q)Z(p+q)}{p^{2}q^{4}(p+q)^{2}}D_{T}^{A\overline{c}c}(-p-q;p,q)(p^{2}\,q^{2}-(p\cdot
q)^{2}))C^{AAA}(-p,p+q,-q),$ (25b) $\displaystyle K^{(c)}$
$\displaystyle=-\frac{N_{c}g^{2}}{2}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{G(q)G(p+q)}{q^{2}(p+q)^{2}}D_{L}^{A\overline{c}c}(q;-p-q,p)\Gamma^{c^{*}cc}(-p,p+q,-q),$
(25c) $\displaystyle K^{(d)}$
$\displaystyle=-\frac{N_{c}g^{2}}{3}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{G(q)G(p+q)}{p^{2}q^{4}(p+q)^{2}}(p^{2}\,q^{2}-(p\cdot
q)^{2}))D_{T}^{A\overline{c}c}(-p;p+q,q)\Gamma^{AA^{*}c}(-q,-p,p+q),$ (25d)
$\displaystyle K^{(e_{1})}$
$\displaystyle=\frac{N_{c}g^{2}}{3}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{Z(p+q)G(q)}{p^{2}q^{4}(p+q)^{4}}(q^{2}+2\,p\cdot
q)(3\,p^{4}+(p\cdot q)^{2}+2\,p^{2}(q^{2}+3\,p\cdot q))$
$\displaystyle\quad\times C^{AAA}(-p,p+q,-q)\Gamma^{AA^{*}c}(p+q,-p,-q),$
(25e) $\displaystyle K^{(e_{2})}$
$\displaystyle=\frac{N_{c}g^{2}}{3}\xi\int\frac{d^{4}q}{(2\pi)^{4}}\frac{G(q)}{q^{4}(p+q)^{4}}(p^{2}q^{2}-(p\cdot
q)^{2})C^{AAA}(-p,p+q,-q)\Gamma^{AA^{*}c}(p+q,-p,-q).$ (25f)

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_{\mu}q_{\nu}r_{\rho}\Gamma^{AAA}_{\mu\nu\rho}(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
$\xi$.

Figure 12: Contributions of individual diagrams at two values of $\xi$ for the
ghost (left) and the gluon (right) NIs.

##  Appendix B Renormalization group properties of correlation functions

Here we provide the RG rescalings of the correlation functions and couplings
${\cal O}_{i}$ used in the NIs, schematically provided in (16).

Since all quantities in the identities are renormalized ones, they behave
under a change of the renormalization group scale $\mu\rightarrow\nu$ as
follows:

$\displaystyle D_{\mu\nu}(p^{2},\mu^{2})$
$\displaystyle=\frac{Z_{3}(\Lambda,\nu)}{Z_{3}(\Lambda,\mu)}D_{\mu\nu}(p^{2},\nu^{2}),$
$\displaystyle D_{G}(p^{2},\mu^{2})$
$\displaystyle=\frac{\widetilde{Z}_{3}(\Lambda,\nu)}{\widetilde{Z}_{3}(\Lambda,\mu)}D_{G}(p^{2},\nu^{2})\,,$
$\displaystyle D^{AAA}(p_{i}^{2},\mu^{2})$
$\displaystyle=\frac{Z_{1}(\Lambda,\mu)}{Z_{1}(\Lambda,\nu)}D^{AAA}(p_{i}^{2},\nu^{2})\,,$
$\displaystyle D^{A\bar{c}c}(p_{i}^{2},\mu^{2})$
$\displaystyle=\frac{\widetilde{Z}_{1}(\Lambda,\mu)}{\widetilde{Z}_{1}(\Lambda,\nu)}D^{A\bar{c}c}(p_{i}^{2},\nu^{2})\,,$
$\displaystyle D^{c*cc}(p_{i}^{2},\mu^{2})$
$\displaystyle=\frac{\widetilde{Z}_{1}(\Lambda,\mu)}{\widetilde{Z}_{1}(\Lambda,\nu)}D^{c*cc}(p_{i}^{2},\nu^{2})\,,$
$\displaystyle D^{A*Ac}(p_{i}^{2},\mu^{2})$
$\displaystyle=\frac{\widetilde{Z}_{1}(\Lambda,\mu)}{\widetilde{Z}_{1}(\Lambda,\nu)}D^{A*Ac}(p_{i}^{2},\nu^{2})\,,$
$\displaystyle g(\mu)$
$\displaystyle=\frac{Z_{g}(\Lambda,\nu)}{Z_{g}(\Lambda,\mu)}g(\nu)\,,$
$\displaystyle\xi(\mu)$
$\displaystyle=\frac{Z_{3}(\Lambda,\nu)}{Z_{3}(\Lambda,\mu)}\xi(\nu)\,,$ (26)

where $Z_{3}$, $\widetilde{Z}_{3}$, $Z_{1}$, $\widetilde{Z}_{1}$, and $Z_{g}$
are the renormalization constants for the gluon propagator, the ghost
propagator, the three-gluon vertex, the ghost-gluon vertex, and the coupling
which are related by the STIs

$\displaystyle
Z_{1}^{2}=Z_{g}^{2}Z_{3}^{3},\quad\widetilde{Z}_{1}^{2}=Z_{g}^{2}Z_{3}\widetilde{Z}_{3}^{2}\,.$
(27)

##  Appendix C Consistent UV limit of correlation functions

Here we discuss the self-consistency of the UV behavior of the NIs. The
employed approximation is exact at the perturbative one-loop level.
Consequently, on the right hand side the correct $\xi$-dependent part of the
anomalous dimension must emerge. This can be seen as follows. Consider the
loop integrals with bare dressing functions for large loop momenta $q$. The
angle integrals can then be performed and lead to

$\displaystyle(a)$
$\displaystyle\xrightarrow{UV}\frac{3\omega}{88}\ln\left(\frac{p^{2}}{\mu^{2}}\right)\,,$
$\displaystyle(b)$
$\displaystyle\xrightarrow{UV}-\frac{9\omega}{88}\ln\left(\frac{p^{2}}{\mu^{2}}\right)\,,$
$\displaystyle(c)$
$\displaystyle\xrightarrow{UV}\frac{3\omega}{22}\ln\left(\frac{p^{2}}{\mu^{2}}\right)\,,$
$\displaystyle(d)$
$\displaystyle\xrightarrow{UV}-\frac{3\omega}{44}\ln\left(\frac{p^{2}}{\mu^{2}}\right)\,,$
$\displaystyle(e_{1})$
$\displaystyle\xrightarrow{UV}\frac{9\omega}{44}\ln\left(\frac{p^{2}}{\mu^{2}}\right)\,,$
$\displaystyle(e_{2})$ $\displaystyle\xrightarrow{UV}0\,.$ (28)

Only the logarithmic parts were kept and $\mu$ is a renormalization scale.
Summing up the corresponding coefficients leads to

$\displaystyle\partial_{\xi}\delta=\frac{3}{44}\,,\qquad\partial_{\xi}\gamma=\frac{3}{22}\,,$
(29)

as can be checked with Tab. 1. This is consistent to one-loop order with the
left-hand side of the equation.

| Anomalous dimension
---|---
Ghost propagator | $\delta=-\frac{9-3\xi}{44}$
Gluon propagator | $\gamma=-\frac{13-3\xi}{22}$
Ghost-gluon vertex | $\gamma^{ghg}=-\frac{3\xi}{22}$
Three-gluon vertex | $\gamma^{\text{3g}}=\frac{17-9\xi}{44}$
Table 1: The one-loop anomalous dimensions of the propagators and vertices.

##  Appendix D Dependence on the parameter $\beta_{0}$

The model we employ for the antifield vertices depends on one parameter
$\beta_{0}$. For the results shown in the main part we used $\beta_{0}=1$. We
tested the influence of this parameter by calculating the propagators also
with $\beta_{0}=0.5$. The model function $H(\bar{p}^{2})$ for these two
choices is shown in Fig. 13. In the quantitatively relevant regime around
$1\,\text{GeV}$ the two parameter values lead only to small differences for
all $\xi$. The rise in the UV for larger $\xi$ comes directly from the
anomalous dimension of the antifield vertices.

The propagators obtained from $\beta_{0}=0.5$ are compared to the ones from
$\beta_{0}=1$ in Fig. 14. We can see that changing $\beta_{0}$ affects
basically only the IR. Only for $\xi>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.

Figure 13: The model employed for the antifield vertices for two values of the
parameter $\beta_{0}$ at different values for $\xi$.

Figure 14: Ghost dressing function (top) and gluon propagator (bottom) for
$\beta_{0}=0.5$ and $1$ at various value of $\xi$.

## References

  * [1] M. Q. Huber, “Gluon and ghost propagators in linear covariant gauges”, Phys. Rev. D91 no. 8, (2015) 085018, arXiv:1502.04057 [hep-ph].
  * [2] M. Q. Huber, “Correlation functions of Landau gauge Yang-Mills theory”, Phys. Rev. D 101 no. 11, (2020) 114009, arXiv:2003.13703 [hep-ph].
  * [3] M. Q. Huber, “Nonperturbative properties of Yang–Mills theories”, Phys. Rept. 879 (2020) 1–92, arXiv:1808.05227 [hep-ph].
  * [4] A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, “Landau gauge Yang-Mills correlation functions”, Phys. Rev. D94 no. 5, (2016) 054005, arXiv:1605.01856 [hep-ph].
  * [5] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, “The nonperturbative functional renormalization group and its applications”, Phys. Rept. 910 (2021) 1–114, arXiv:2006.04853 [cond-mat.stat-mech].
  * [6] A. Cucchieri and T. Mendes, “What’s up with IR gluon and ghost propagators in Landau gauge? A puzzling answer from huge lattices”, PoS LAT2007 (2007) 297, arXiv:0710.0412 [hep-lat].
  * [7] A. Cucchieri and T. Mendes, “Constraints on the IR behavior of the ghost propagator in Yang-Mills theories”, Phys. Rev. D78 (2008) 094503, arXiv:0804.2371 [hep-lat].
  * [8] I. L. Bogolubsky, E. M. Ilgenfritz, M. Müller-Preussker, and A. Sternbeck, “The Landau gauge gluon and ghost propagators in 4D SU(3) gluodynamics in large lattice volumes”, PoS LAT2007 (2007) 290, arXiv:0710.1968 [hep-lat].
  * [9] I. L. Bogolubsky, E. M. Ilgenfritz, M. Müller-Preussker, and A. Sternbeck, “Lattice gluodynamics computation of Landau gauge Green’s functions in the deep infrared”, Phys. Lett. B676 (2009) 69–73, arXiv:0901.0736 [hep-lat].
  * [10] R. Alkofer and L. von Smekal, “The infrared behavior of QCD Green’s functions: Confinement, dynamical symmetry breaking, and hadrons as relativistic bound states”, Phys. Rept. 353 (2001) 281, arXiv:hep-ph/0007355.
  * [11] C. S. Fischer, “Infrared properties of QCD from Dyson-Schwinger equations”, J. Phys. G32 (2006) R253–R291, arXiv:hep-ph/0605173.
  * [12] C. S. Fischer, A. Maas, and J. M. Pawlowski, “On the infrared behavior of Landau gauge Yang-Mills theory”, Annals Phys. 324 (2009) 2408–2437, arXiv:0810.1987 [hep-ph].
  * [13] A. Aguilar, D. Binosi, and J. Papavassiliou, “Gluon and ghost propagators in the Landau gauge: Deriving lattice results from Schwinger-Dyson equations”, Phys.Rev. D78 (2008) 025010, arXiv:0802.1870 [hep-ph].
  * [14] D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel, and H. Verschelde, “A refinement of the Gribov-Zwanziger approach in the Landau gauge: infrared propagators in harmony with the lattice results”, Phys. Rev. D78 (2008) 065047, arXiv:0806.4348 [hep-th].
  * [15] D. Dudal, S. Sorella, and N. Vandersickel, “The dynamical origin of the refinement of the Gribov-Zwanziger theory”, Phys.Rev. D84 (2011) 065039, arXiv:1105.3371 [hep-th].
  * [16] A. Maas, “Describing gauge bosons at zero and finite temperature”, Phys.Rept. 524 (2013) 203–300, arXiv:1106.3942 [hep-ph].
  * [17] P. Boucaud, J. P. Leroy, A. L. Yaouanc, J. Micheli, O. Pene, and J. Rodriguez-Quintero, “The Infrared Behaviour of the Pure Yang-Mills Green Functions”, Few Body Syst. 53 (2012) 387–436, arXiv:1109.1936 [hep-ph].
  * [18] M. Quandt, H. Reinhardt, and J. Heffner, “Covariant variational approach to Yang-Mills theory”, Phys. Rev. D89 no. 6, (2014) 065037, arXiv:1310.5950 [hep-th].
  * [19] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, “Gluon mass scale through nonlinearities and vertex interplay”, Phys. Rev. D100 no. 9, (2019) 094039, arXiv:1909.09826 [hep-ph].
  * [20] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, “A window on infrared QCD with small expansion parameters”, arXiv:2106.04526 [hep-th].
  * [21] A. Aguilar, D. Binosi, and J. Papavassiliou, “Yang-Mills two-point functions in linear covariant gauges”, Phys.Rev. D91 (2015) 085014, arXiv:1501.07150 [hep-ph].
  * [22] J. C. Breckenridge, M. J. Lavelle, and T. G. Steele, “The Nielsen identities for the two point functions of QED and QCD”, Z. Phys. C 65 (1995) 155–164, arXiv:hep-th/9407028.
  * [23] R. Alkofer, C. S. Fischer, H. Reinhardt, and L. von Smekal, “On the infrared behaviour of gluons and ghosts in ghost- antighost symmetric gauges”, Phys. Rev. D68 (2003) 045003, arXiv:hep-th/0304134.
  * [24] A. Aguilar and J. Papavassiliou, “Infrared finite ghost propagator in the Feynman gauge”, Phys.Rev. D77 (2008) 125022, arXiv:0712.0780 [hep-ph].
  * [25] M. Q. Huber, K. Schwenzer, and R. Alkofer, “On the infrared scaling solution of SU(N) Yang-Mills theories in the maximally Abelian gauge”, Eur. Phys. J. C68 (2010) 581–600, arXiv:0904.1873 [hep-th].
  * [26] A. C. Aguilar, D. Binosi, and J. Papavassiliou, “Schwinger mechanism in linear covariant gauges”, Phys. Rev. D 95 no. 3, (2017) 034017, arXiv:1611.02096 [hep-ph].
  * [27] T. De Meerleer, D. Dudal, S. P. Sorella, P. Dall’Olio, and A. Bashir, “Landau-Khalatnikov-Fradkin Transformations, Nielsen Identities, Their Equivalence and Implications for QCD”, Phys. Rev. D 101 no. 8, (2020) 085005, arXiv:1911.01907 [hep-th].
  * [28] R. F. Sobreiro and S. P. Sorella, “A study of the Gribov copies in linear covariant gauges in Euclidean Yang-Mills theories”, JHEP 06 (2005) 054, arXiv:hep-th/0506165.
  * [29] M. A. L. Capri, D. Dudal, D. Fiorentini, M. S. Guimaraes, I. F. Justo, A. D. Pereira, B. W. Mintz, L. F. Palhares, R. F. Sobreiro, and S. P. Sorella, “Exact nilpotent nonperturbative BRST symmetry for the Gribov-Zwanziger action in the linear covariant gauge”, Phys. Rev. D 92 no. 4, (2015) 045039, arXiv:1506.06995 [hep-th].
  * [30] M. A. L. Capri, D. Fiorentini, M. S. Guimaraes, B. W. Mintz, L. F. Palhares, S. P. Sorella, D. Dudal, I. F. Justo, A. D. Pereira, and R. F. Sobreiro, “More on the nonperturbative Gribov-Zwanziger quantization of linear covariant gauges”, Phys. Rev. D 93 no. 6, (2016) 065019, arXiv:1512.05833 [hep-th].
  * [31] M. A. L. Capri, D. Dudal, A. D. Pereira, D. Fiorentini, M. S. Guimaraes, B. W. Mintz, L. F. Palhares, and S. P. Sorella, “Nonperturbative aspects of Euclidean Yang-Mills theories in linear covariant gauges: Nielsen identities and a BRST-invariant two-point correlation function”, Phys. Rev. D 95 no. 4, (2017) 045011, arXiv:1611.10077 [hep-th].
  * [32] F. Siringo, “Gluon propagator in Feynman gauge by the method of stationary variance”, Phys.Rev. D90 no. 9, (2014) 094021, arXiv:1408.5313 [hep-ph].
  * [33] F. Siringo, “Second order gluon polarization for SU(N) theory in a linear covariant gauge”, Phys. Rev. D92 no. 7, (2015) 074034, arXiv:1507.00122 [hep-ph].
  * [34] F. Siringo and G. Comitini, “The gluon propagator in linear covariant $R_{\xi}$ gauges”, arXiv:1806.08397 [hep-ph].
  * [35] A. Cucchieri, T. Mendes, and E. M. Santos, “Covariant gauge on the lattice: A New implementation”, Phys.Rev.Lett. 103 (2009) 141602, arXiv:0907.4138 [hep-lat].
  * [36] A. Cucchieri, T. Mendes, G. M. Nakamura, and E. M. S. Santos, “Gluon Propagators in Linear Covariant Gauge”, PoS FACESQCD (2010) 026, arXiv:1102.5233 [hep-lat].
  * [37] P. Bicudo, D. Binosi, N. Cardoso, O. Oliveira, and P. J. Silva, “Lattice gluon propagator in renormalizable $\xi$ gauges”, Phys. Rev. D92 no. 11, (2015) 114514, arXiv:1505.05897 [hep-lat].
  * [38] A. Cucchieri, D. Dudal, T. Mendes, O. Oliveira, M. Roelfs, and P. J. Silva, “Lattice Computation of the Ghost Propagator in Linear Covariant Gauges”, PoS LATTICE2018 (2018) 252, arXiv:1811.11521 [hep-lat].
  * [39] A. Cucchieri, D. Dudal, T. Mendes, O. Oliveira, M. Roelfs, and P. J. Silva, “Faddeev-Popov Matrix in Linear Covariant Gauge: First Results”, Phys. Rev. D 98 no. 9, (2018) 091504, arXiv:1809.08224 [hep-lat].
  * [40] A. Cucchieri, A. Maas, and T. Mendes, “Infrared-suppressed gluon propagator in 4-D Yang-Mills theory in a Landau-like gauge”, Mod. Phys. Lett. A22 (2007) 2429–2438, arXiv:hep-lat/0701011 [hep-lat].
  * [41] A. Andrasi and J. C. Taylor, “Renormalization in a Landau to Coulomb interpolating gauge in Yang-Mills theory”, arXiv:2103.01683 [hep-th].
  * [42] D. Dudal, J. Gracey, V. Lemes, M. Sarandy, R. F. Sobreiro, S. P. Sorella, and H. Verschelde, “An analytic study of the off-diagonal mass generation for Yang-Mills theories in the maximal Abelian gauge”, Phys. Rev. D70 (2004) 114038, arXiv:hep-th/0406132.
  * [43] M. Q. Huber, On gauge fixing aspects of the infrared behavior of Yang-Mills Green functions.  Springer, 2010.  arXiv:1005.1775 [hep-th].  Ph.D.Thesis, Karl-Franzens-University Graz, Theses series.
  * [44] M. A. L. Capri, R. F. Sobreiro, and S. P. Sorella, “Interpolating among the Landau, Coulomb and maximal Abelian gauges”, Phys. Rev. D73 (2006) 041701, arXiv:hep-th/0512096.
  * [45] M. A. L. Capri, R. F. Sobreiro, S. P. Sorella, and R. Thibes, “Renormalizability of a generalized gauge fixing interpolating among the Coulomb, Landau and maximal Abelian gauges”, Annals Phys. 322 (2007) 1776–1789, arXiv:hep-th/0607117 [hep-th].
  * [46] D. Dudal, J. A. Gracey, V. E. R. Lemes, R. F. Sobreiro, S. P. Sorella, R. Thibes, and H. Verschelde, “Remarks on a class of renormalizable interpolating gauges”, JHEP 07 (2005) 059, arXiv:hep-th/0505037 [hep-th].
  * [47] C. Becchi, A. Rouet, and R. Stora, “Renormalization of Gauge Theories”, Annals Phys. 98 (1976) 287–321.
  * [48] I. V. Tyutin, “Gauge Invariance in Field Theory and Statistical Physics in Operator Formalism”, LEBEDEV-75-39 (1975) , arXiv:0812.0580 [hep-th].
  * [49] I. A. Batalin and G. A. Vilkovisky, “Gauge Algebra and Quantization”, Phys. Lett. B 102 (1981) 27–31.
  * [50] I. A. Batalin and G. a. Vilkovisky, “FEYNMAN RULES FOR REDUCIBLE GAUGE THEORIES”, Phys. Lett. B 120 (1983) 166–170.
  * [51] I. A. Batalin and G. A. Vilkovisky, “Quantization of Gauge Theories with Linearly Dependent Generators”, Phys. Rev. D 28 (1983) 2567–2582.
  * [52] A. Fuster, M. Henneaux, and A. Maas, “BRST quantization: A Short review”, Int. J. Geom. Meth. Mod. Phys. 2 (2005) 939–964, arXiv:hep-th/0506098.
  * [53] G. Eichmann, R. Williams, R. Alkofer, and M. Vujinovic, “The three-gluon vertex in Landau gauge”, Phys.Rev. D89 (2014) 105014, arXiv:1402.1365 [hep-ph].
  * [54] D. Binosi and J. Papavassiliou, “Pinch technique and the Batalin-Vilkovisky formalism”, Phys. Rev. D 66 (2002) 025024, arXiv:hep-ph/0204128.
  * [55] D. Binosi, J. Papavassiliou, and A. Pilaftsis, “Displacement operator formalism for renormalization and gauge dependence to all orders”, Phys. Rev. D 71 (2005) 085007, arXiv:hep-ph/0501259.
  * [56] A. Sternbeck, arXiv:hep-lat/0609016, PhD thesis, Humboldt-Universität zu Berlin, 2006\.
  * [57] A. Maas, “Constraining the gauge-fixed Lagrangian in minimal Landau gauge”, arXiv:1907.10435 [hep-lat].
  * [58] A. Cucchieri, A. Maas, and T. Mendes, “Three-point vertices in Landau-gauge Yang-Mills theory”, Phys. Rev. D77 (2008) 094510, arXiv:0803.1798 [hep-lat].
  * [59] A. Sternbeck, P.-H. Balduf, A. Kızılersu, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. G. Williams, “Triple-gluon and quark-gluon vertex from lattice QCD in Landau gauge”, PoS LATTICE2016 (2017) 349, arXiv:1702.00612 [hep-lat].
  * [60] A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, “On the zero crossing of the three-gluon vertex”, Phys. Lett. B761 (2016) 444–449, arXiv:1607.01278 [hep-ph].
  * [61] A. G. Duarte, O. Oliveira, and P. J. Silva, “Further Evidence For Zero Crossing On The Three Gluon Vertex”, Phys. Rev. D94 no. 7, (2016) 074502, arXiv:1607.03831 [hep-lat].
  * [62] A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, and J. Rodríguez-Quintero, “Infrared facets of the three-gluon vertex”, Phys. Lett. B 818 (2021) 136352, arXiv:2102.04959 [hep-ph].
  * [63] A. Blum, M. Q. Huber, M. Mitter, and L. von Smekal, “Gluonic three-point correlations in pure Landau gauge QCD”, Phys. Rev. D 89 (2014) 061703(R), arXiv:1401.0713 [hep-ph].
  * [64] A. Aguilar, D. Binosi, D. Ibáñez, and J. Papavassiliou, “Effects of divergent ghost loops on the Green’s functions of QCD”, Phys.Rev. D89 (2014) 085008, arXiv:1312.1212 [hep-ph].
  * [65] M. Q. Huber.  https://github.com/markusqh/YM_data. Github repository.
  * [66] L. von Smekal, A. Hauck, and R. Alkofer, “A solution to coupled Dyson-Schwinger equations for gluons and ghosts in Landau gauge”, Ann. Phys. 267 (1998) 1, arXiv:hep-ph/9707327.
  * [67] L. von Smekal, R. Alkofer, and A. Hauck, “The infrared behavior of gluon and ghost propagators in Landau gauge QCD”, Phys. Rev. Lett. 79 (1997) 3591–3594, arXiv:hep-ph/9705242.
  * [68] D. Zwanziger, “Nonperturbative Landau gauge and infrared critical exponents in QCD”, Phys.Rev. D65 (2002) 094039, arXiv:hep-th/0109224 [hep-th].
  * [69] C. Lerche and L. von Smekal, “On the infrared exponent for gluon and ghost propagation in Landau gauge QCD”, Phys.Rev. D65 (2002) 125006, arXiv:hep-ph/0202194.
  * [70] J. M. Pawlowski, D. F. Litim, S. Nedelko, and L. von Smekal, “Infrared behaviour and fixed points in Landau gauge QCD”, Phys. Rev. Lett. 93 (2004) 152002, arXiv:hep-th/0312324.
  * [71] R. Alkofer, M. Q. Huber, and K. Schwenzer, “Infrared singularities in Landau gauge Yang-Mills theory”, Phys. Rev. D81 (2010) 105010, arXiv:0801.2762 [hep-th].
  * [72] M. Q. Huber and M. Mitter, “CrasyDSE: A Framework for solving Dyson-Schwinger equations”, Comput.Phys.Commun. 183 (2012) 2441–2457, arXiv:1112.5622 [hep-th].
  * [73] M. Q. Huber and L. von Smekal, “Spurious divergences in Dyson-Schwinger equations”, JHEP 1406 (2014) 015, arXiv:1404.3642 [hep-ph].
  * [74] A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, “Unified description of seagull cancellations and infrared finiteness of gluon propagators”, Phys. Rev. D94 no. 4, (2016) 045002, arXiv:1604.08456 [hep-ph].
  * [75] A. I. Davydychev, P. Osland, and O. V. Tarasov, “Two loop QCD vertices, WST identities and RG quantities” in Workshop on the Exact Renormalization Group.  11, 1998.  arXiv:hep-ph/9811503.
  * [76] A. I. Davydychev, P. Osland, and O. V. Tarasov, “Two loop three gluon vertex in zero momentum limit”, Phys. Rev. D 58 (1998) 036007, arXiv:hep-ph/9801380.
  * [77] K. Chetyrkin and A. Retey, “Three loop three linear vertices and four loop similar to MOM beta functions in massless QCD”, arXiv:hep-ph/0007088.
  * [78] A. C. Aguilar, D. Binosi, J. Papavassiliou, and J. Rodriguez-Quintero, “Non-perturbative comparison of QCD effective charges”, Phys. Rev. D 80 (2009) 085018, arXiv:0906.2633 [hep-ph].
  * [79] P. Boucaud et al., “IR finiteness of the ghost dressing function from numerical resolution of the ghost SD equation”, JHEP 06 (2008) 012, arXiv:0801.2721 [hep-ph].
  * [80] A. Maas, “Constructing non-perturbative gauges using correlation functions”, Phys. Lett. B689 (2010) 107–111, arXiv:0907.5185 [hep-lat].
  * [81] A. Sternbeck and M. Müller-Preussker, “Lattice evidence for the family of decoupling solutions of Landau gauge Yang-Mills theory”, Phys.Lett. B726 (2013) 396–403, arXiv:1211.3057 [hep-lat].
  * [82] M. Q. Huber, C. S. Fischer, and H. Sanchis-Alepuz, “Spectrum of scalar and pseudoscalar glueballs from functional methods”, Eur. Phys. J. C 80 no. 11, (2020) 1077, arXiv:2004.00415 [hep-ph].
  * [83] W. Kern, M. Q. Huber, and R. Alkofer, “Spectral dimension as a tool for analyzing nonperturbative propagators”, Phys. Rev. D 100 no. 9, (2019) 094037, arXiv:1909.12207 [hep-th].
  * [84] M. Q. Huber, R. Alkofer, and K. Schwenzer, “Infrared scaling solutions beyond the Landau gauge: The maximally Abelian gauge and Abelian infrared dominance”, PoS FACESQCD (2010) 001, arXiv:1103.0236 [hep-th].