# Gauge dependence of the quark gap equation: an exploratory study

José Roberto Lessa  Laboratório de Física Teórica e Computacional,
Universidade Cidade de São Paulo, Rua Galvão Bueno 868, São Paulo, São Paulo
01506-000, Brazil<EMAIL_ADDRESS>Fernando E. Serna
<EMAIL_ADDRESS>Laboratório de Física Teórica e
Computacional, Universidade Cidade de São Paulo, Rua Galvão Bueno 868, São
Paulo, São Paulo 01506-000, Brazil  Departamento de Física, Universidad de
Sucre, Carrera 28 No. 5-267, Barrio Puerta Roja, Sincelejo 700001, Colombia
Bruno El-Bennich<EMAIL_ADDRESS>Laboratório de Física
Teórica e Computacional, Universidade Cidade de São Paulo, Rua Galvão Bueno
868, São Paulo, São Paulo 01506-000, Brazil  Departamento de Física, ICAFQ,
Universidade Federal de São Paulo, Diadema, São Paulo 09913-030, Brazil  Adnan
Bashir  adnan.bashir<EMAIL_ADDRESS>Instituto de Física y
Matemáticas, Universidad Michoacana de San Nicolás de Hidalgo, Morelia,
Michoacán 58040, Mexico  Theory Center, Jefferson Lab, Newport News, VA 23606,
USA  Orlando Oliveira<EMAIL_ADDRESS>CFisUC, Department of Physics,
University of Coimbra, 3004 516 Coimbra, Portugal

###### Abstract

We study the gauge dependence of the quark propagator in quantum
chromodynamics (QCD) by solving the gap equation with a nonperturbative quark-
gluon vertex which is constrained by longitudinal and transverse Slavnov-
Taylor identities, the discrete charge conjugation and parity symmetries and
which is free of kinematic singularities in the limit of equal incoming and
outgoing quark momenta. We employ gluon propagators in renormalizable
$R_{\xi}$ gauges obtained in lattice QCD studies. We report the dependence of
the nonperturbative quark propagator on the gauge parameter, in particular we
observe an increase, proportional to the gauge-fixing parameter, of the mass
function in the infrared domain, whereas the wave renormalization decreases
within the range $0\leq\xi\leq 1$ considered here. The chiral quark condensate
reveals a mild gauge dependence in the region of $\xi$ investigated. We
comment on how to build and improve upon this exploratory study in future in
conjunction with generalized gauge covariance relations for QCD.

##  I Introduction

Gauge symmetry and its innumerable consequences have played a fundamental role
in the development of modern quantum field theory during the past century.
Amongst its celebrated implications, the Ward-Takahashi identities (WTI) Ward
(1950); Takahashi (1957) in quantum electrodynamics (QED) and the
corresponding Slavnov-Taylor identities (STI) Slavnov (1972); Taylor (1971) of
quantum chromodynamics (QCD) relate different $n$-point Green functions to
each other. In particular, the WTI and the STI connecting the divergence of
the fermion-boson vertex to the fermion propagator help us identify the
longitudinal part Ball and Chiu (1980); Aguilar et al. (2017a); Bermudez et
al. (2017) of this three-point vertex, whether it be the quark-photon or the
quark-gluon vertex. These are exact nonperturbative relations which are
observed order by order in perturbation theory. There exists a plethora of
work within the nonperturbative exploration of Dyson-Schwinger equations
(DSEs) which incorporates these identities in studying dynamical chiral
symmetry breaking (DCSB) via the electron/quark gap equation.

Whereas the usual WTI or STI relates the divergence of the three-point
fermion-boson vertex to the inverse fermion propagator, there exist transverse
Takahashi identities (TTI) and transverse Slavnov-Taylor identities (TSTI)
which play a similar role for the curl of the fermion-boson vertex Takahashi
(1985); Kondo (1997); He et al. (2000); He (2006, 2007). The TTI and TSTI are
richer and more complicated in their structure, and they shed light on the
transverse part of the fermion-boson vertex. They have been employed to
compute the critical coupling and study its gauge independence, as well as the
quark condensate and pion decay constant He (2009); Qin et al. (2013, 2014);
Albino et al. (2019, 2021).

While the longitudinal and transverse gauge identities relate different
$n$-point Green functions with each other, another important consequence of
gauge covariance in QED are the Landau-Khalatnikov-Fradkin transformations
(LKFTs) which describe how the individual Green functions respond to an
arbitrary gauge transformation Landau and Khalatnikov (1955); Fradkin (1955).
The LKFTs are a well-defined set of transformations which leave the DSEs and
related WTIs of the fermion-boson vertex form invariant. While the STIs are
the QCD generalization of the WTIs, the equivalent generalization of the LKFTs
has only recently been derived by two groups employing two different methods:
(i) ABG Aslam et al. (2016) through direct generalization of the method
employed by Landau and Khalatnikov and (ii) MDSDB De Meerleer et al. (2018,
2020) based on the introduction of a gauge invariant transverse gauge field.
Moreover, the Nielsen identities allow to study the variation of Green
functions as derivatives with respect to the gauge parameter De Meerleer et
al. (2020); Nielsen (1975). Nevertheless, even with these formal results
available, starting from the knowledge of Green functions in a given gauge,
their explicit extraction in another gauge remains a nontrivial problem in QCD
even at one loop in perturbation theory. This is already the case for two-
point Green functions, let alone their nonperturbative transformation and that
of three- and higher $n$-point Green functions.

These local transformations are in no simple manner amenable to
straightforward comparisons of Green functions in different gauges, especially
in momentum space, nor to explicitly prove the gauge invariance of physical
observables. However, it has explicitly been shown in QED, and verified in
QED3 by constructing a nonperturbative vertex Bashir et al. (2011), studying
gauge-covariance relations Bashir et al. (1999, 2000) and generating numerical
solutions of the electron gap equation, that the electron condensate is
manifestly gauge independent once its propagator is correctly LKF transformed
from a given covariant gauge to other covariant gauges Bashir and Raya (2005,
2007). Moreover, the ABG and the MDSDB generalizations of LKFT in QCD confirm
the formal gauge-invariance of the quark condensate after the quark propagator
has been adequately gauge transformed.

We are not able to directly verify this result by gauge transforming every
Green function computed in Landau gauge to repeat the calculation in any other
gauge—it remains a prohibitively difficult task for the time being. We thus
follow a more modest approach to study the gauge dependence of the quark
propagator in QCD and _solve_ the associated DSE, employing quenched gluon
propagators in $R_{\xi}$ gauges from lattice QCD Bicudo et al. (2015), and the
dressed quark-gluon vertex constructed in Ref. Albino et al. (2021), invoking
the STI and TSTI. Our approach not only allows for a comparison of the DSE
solutions in $R_{\xi}$ gauges, i.e., the mass and wave renormalization
functions of the quark for gauge parameters in the range $\xi\in[0,0.5]$,
extrapolated to Feynman gauge, $\xi=1$, but also to compute the quark
condensate as a function of the gauge parameter within this interval. We
emphasize that this is an initial, exploratory study which we expect to shed
light on how far we are from obtaining the gauge independence of the
condensate. This invariance would be obtained through formal incorporation of
the local gauge transformations of every Green function in the gap equation,
i.e., the gluon propagator, the quark-gluon vertex as well as the ghost
propagator and the ghost-quark scattering kernel.

Within the limitations of this hybrid approach, we demonstrate how the
numerical solutions of the quark propagator vary as function of the gauge
parameter in $R_{\xi}$ gauges, and that this behavior leads to a slightly
gauge-dependent quark condensate up to Feynman gauge within the error
estimates involved. Nonetheless, to our knowledge, we present the first DSE
solutions for the quark obtained with all 12 vector structures of the
nonperturbative quark-gluon vertex and the gluon and ghost propagators in
covariant $R_{\xi}$ gauges. Our findings are encouraging and represent an
important first step as they incorporate the gauge covariance into the DSEs
through the formal implementation of gauge identities for the Green functions.
This is a necessary, though not sufficient requirement for calculating gauge-
independent bound-state properties.

This article is organized as follows: in Sec. II we describe our functional
approach to QCD and discuss the quark-gluon vertex, the lattice-extracted
gluon and ghost propagators and the quark-ghost form factor in different
gauges. In Sec. III, we solve the quark gap equation and obtain the mass and
wave renormalization functions in different gauges and calculate the quark
condensate. Final remarks are given in Sect. IV.

##  II Dyson-Schwinger equation in Rξ gauges

The DSEs are the relativistic equations of motion in quantum field theory, see
e.g. Ref. Bashir et al. (2012a) for a review. For a given flavor and in
$R_{\xi}$ gauge, the DSE of the inverse quark propagator in Euclidean space
reads,

$\displaystyle S^{-1}_{\xi}(p)$ $\displaystyle=\,Z_{2}\,i\,\gamma\cdot
p+Z_{4}\,m(\mu)$ $\displaystyle+\,Z_{1}$ $\displaystyle
4\pi\alpha_{s}^{\xi}\\!\int^{\Lambda}\\!\\!\frac{d^{4}k}{(2\pi)^{4}}\
\Delta^{ab}_{\mu\nu}(q)\,\gamma_{\mu}t^{a}\,S_{\xi}(k)\,\Gamma_{\nu}^{b\xi}(k,p)\,,$
(1)

where $m(\mu)$ is the renormalized current-quark mass, $Z_{1}(\mu,\Lambda)$,
$Z_{2}(\mu,\Lambda)$ and $Z_{4}(\mu,\Lambda)$ are the vertex, wave function
and mass renormalization constants, respectively, while
$\Gamma^{a\xi}_{\mu}(k,p)=\Gamma_{\mu}^{\xi}(k,p)t^{a}$ is the dressed quark-
gluon vertex and $t^{a}=\lambda^{a}/2$ are the SU(3) group generators in the
fundamental representation. The gluon propagator in $R_{\xi}$ gauge with
momentum $q=k-p$,

$\Delta^{ab}_{\mu\nu}(q)=\delta^{ab}\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}}\right)\\!\Delta_{\xi}(q^{2})+\delta^{ab}\,\xi\dfrac{q_{\mu}q_{\nu}}{q^{4}}\
$ (2)

is characterized by a nonperturbative transverse dressing function,
$\Delta_{\xi}(\mu^{2})=1/\mu^{2}$, and was studied with different lattice and
functional approaches, e.g. Refs. Bicudo et al. (2015); Huber (2015);
Napetschnig et al. (2021).

The most general Poincaré-covariant form of the solutions to Eq. (II) is
written in terms of covariant scalar and vector amplitudes:

$S_{\xi}(p)=\frac{1}{i\gamma\cdot
p\,A_{\xi}(p^{2})+B_{\xi}(p^{2})}=\,\frac{Z_{\xi}(p^{2})}{i\gamma\cdot
p+M_{\xi}(p^{2})}\ .$ (3)

In the self-energy integral, $\Lambda$ is a Poincaré-invariant ultraviolet
cutoff and $\mu$ is the renormalization scale imposed such that
$\Lambda\gg\mu$. This scale is implicit in our convenient notation:
$A(p^{2})\equiv A(p^{2},\mu^{2})$ and $B(p^{2})\equiv B(p^{2},\mu^{2})$, as is
a flavor index $f$ for these quantities and for all renormalization constants.
The flavor- and gauge-dependent mass and wave renormalization functions are,
respectively:

$\displaystyle M_{\xi}(p^{2})$
$\displaystyle=B_{\xi}(p^{2},\mu^{2})/A_{\xi}(p^{2},\mu^{2})\ ,$ (4)
$\displaystyle Z_{\xi}(p^{2},\mu^{2})$
$\displaystyle=1/A_{\xi}(p^{2},\mu^{2})\ .$ (5)

We choose the renormalization scale, $\mu=4.3$ GeV, for a twofold reason: _i_)
we consistently renormalize the DSE at the scale at which the transverse
dressing function $\Delta_{\xi}(q^{2})$ in $R_{\xi}$ gauge is renormalized
Bicudo et al. (2015) and _ii_) we can compare the dressed functions
$M_{\xi}(p^{2})$ and $Z_{\xi}(p^{2},\Lambda^{2})$ with the solutions of
lattice regularized QCD at this scale Bowman et al. (2004); Oliveira et al.
(2019). We therefore impose Rojas et al. (2013); Serna et al. (2019):
$Z_{\xi}(\mu^{2})=1$ and $M_{\xi}(\mu^{2})\equiv m(\mu^{2})=25$ MeV.

For the dressed quark-gluon vertex we employ the decomposition detailed, e.g.,
in Ref. Albino et al. (2019),

$\displaystyle\Gamma_{\mu}^{\xi}(k,p)$
$\displaystyle=\Gamma^{L\xi}_{\mu}(k,p)+\Gamma^{T\xi}_{\mu}(k,p)$
$\displaystyle=\sum_{i=1}^{4}$
$\displaystyle\,\lambda_{i}^{\xi}(k,p)L_{\mu}^{i}(k,p)+\sum_{i=1}^{8}\tau_{i}^{\xi}(k,p)T_{\mu}^{i}(k,p)\
,$ (6)

where the transverse vertex $\Gamma^{T\xi}_{\mu}(k,p)$ in Eq. (6) is naturally
defined by $i\,q\cdot\Gamma^{T\xi}(k,p)=0$ . The usual STI Slavnov (1972);
Taylor (1971) constrains the longitudinal vertex, $\Gamma^{L\xi}_{\mu}(k,p)$,
to four independent structures but leaves the transverse components
undetermined. The remaining eight tensor structures can be explored with TSTI
derived from the symmetry transformation that involves the Lorentz
transformation acting on the usual infinitesimal gauge transformation.

The TSTI He (2009) constrain but also couple the vector and axial-vector
vertices. However, as shown in Refs. Qin et al. (2013); Albino et al. (2019),
these can be decoupled and merely one identity for the vector vertex is
sufficient to obtain analytic expressions for the transverse form factors
$\tau_{i}^{\xi}(k,p)$. The full vertex is thus described by the form factors
derived in Refs. Albino et al. (2021); El-Bennich et al. (2022),

$\displaystyle\lambda_{1}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}\,G(q^{2})X_{0}^{\xi}(q^{2})\left[A_{\xi}(k^{2})+A_{\xi}(p^{2})\right]\
,$ (7) $\displaystyle\lambda_{2}^{\xi}(k,p)$ $\displaystyle=$ $\displaystyle
G(q^{2})X_{0}^{\xi}(q^{2})\,\frac{A_{\xi}(k^{2})-A_{\xi}(p^{2})}{k^{2}-p^{2}}\
,$ (8) $\displaystyle\lambda_{3}^{\xi}(k,p)$ $\displaystyle=$ $\displaystyle
G(q^{2})X_{0}^{\xi}(q^{2})\,\frac{B_{\xi}(k^{2})-B_{\xi}(p^{2})}{k^{2}-p^{2}}\
,$ (9) $\displaystyle\lambda_{4}^{\xi}(k,p)$ $\displaystyle=$ $\displaystyle\
0\ ,$ (10)

for the longitudinal part and by

$\displaystyle\tau_{1}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle-\frac{Y_{1}}{2(k^{2}-p^{2})\nabla(k,p)}\ ,$ (11)
$\displaystyle\tau_{2}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle-\frac{Y_{5}-3Y_{3}}{4(k^{2}-p^{2})\nabla(k,p)}\ ,$ (12)
$\displaystyle\tau_{3}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,G(q^{2})X_{0}^{\xi}(q^{2})\left[\frac{A_{\xi}(k^{2})-A_{\xi}(p^{2})}{k^{2}-p^{2}}\right]$
(13) $\displaystyle+$
$\displaystyle\frac{Y_{2}}{4\nabla(k,p)}-\frac{(k+p)^{2}(Y_{3}-Y_{5})}{8(k^{2}-p^{2})\nabla(k,p)}\
,$ $\displaystyle\tau_{4}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle-\frac{6Y_{4}+Y_{6}^{A}}{8\nabla(k,p)}-\frac{(k+p)^{2}Y_{7}^{S}}{8(k^{2}-p^{2})\nabla(k,p)}\
,$ (14) $\displaystyle\tau_{5}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle-G(q^{2})X_{0}^{\xi}(q^{2})\left[\frac{B_{\xi}(k^{2})-B_{\xi}(p^{2})}{k^{2}-p^{2}}\right]$
(15) $\displaystyle-$ $\displaystyle\frac{2Y_{4}+Y_{6}^{A}}{2(k^{2}-p^{2})}\
,$ $\displaystyle\tau_{6}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle\frac{(k-p)^{2}Y_{2}}{4(k^{2}-p^{2})\nabla(k,p)}-\frac{Y_{3}-Y_{5}}{8\nabla(k,p)}\
,$ (16) $\displaystyle\tau_{7}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle\frac{q^{2}(6Y_{4}+Y_{6}^{A})}{4(k^{2}-p^{2})\nabla(k,p)}+\frac{Y_{7}^{S}}{4\nabla(k,p)}\
,$ (17) $\displaystyle\tau_{8}^{\xi}(k,p)$ $\displaystyle=$
$\displaystyle-G(q^{2})X_{0}^{\xi}(q^{2})\left[\frac{A_{\xi}(k^{2})-A_{\xi}(p^{2})}{k^{2}-p^{2}}\right]$
(18) $\displaystyle-$ $\displaystyle\frac{2Y_{8}^{A}}{k^{2}-p^{2}}\ ,$

for the transverse components, where the Gram determinant is defined as:
$\nabla(k,p)=k^{2}p^{2}-(k\cdot p)^{2}$.

$\xi$ |  $Z$ |  $M_{0}^{2}$ [GeV2] |  $M_{1}^{2}$ [GeV2] |  $M_{2}^{2}$ [GeV2] |  $M_{3}^{4}$ [GeV4] |  $\chi^{2}$/dof
---|---|---|---|---|---|---
0.0 |  $8.407\pm 0.031$ |  $0.071\pm 0.023$ |  $2.783\pm 0.128$ |  $0.532\pm 0.055$ |  $0.380\pm 0.012$ |  1.209
0.1 |  $8.321\pm 0.023$ |  $0.071\pm 0.018$ |  $2.712\pm 0.099$ |  $0.517\pm 0.043$ |  $0.374\pm 0.009$ |  0.734
0.2 |  $8.238\pm 0.018$ |  $0.074\pm 0.015$ |  $2.677\pm 0.076$ |  $0.515\pm 0.033$ |  $0.371\pm 0.007$ |  0.421
0.3 |  $8.160\pm 0.012$ |  $0.069\pm 0.010$ |  $2.617\pm 0.051$ |  $0.499\pm 0.023$ |  $0.367\pm 0.005$ |  0.205
0.4 |  $8.072\pm 0.011$ |  $0.076\pm 0.010$ |  $2.595\pm 0.046$ |  $0.499\pm 0.020$ |  $0.363\pm 0.004$ |  0.178
0.5 |  $7.996\pm 0.009$ |  $0.075\pm 0.009$ |  $2.570\pm 0.041$ |  $0.500\pm 0.018$ |  $0.362\pm 0.004$ |  0.162
0.7 |  $7.831$ |  $0.077$ |  $2.488$ |  $0.487$ |  $0.354$ |  —
1.0 |  $7.585$ |  $0.080$ |  $2.378$ |  $0.473$ |  $0.344$ |  —
Table 1: Parameters of the gluon-dressing parametrization (23) as a function
of the gauge parameter $\xi$. Note that the bare gluon propagators of Ref.
Bicudo et al. (2015) were fitted which must be renormalized as
$\Delta_{\xi}(\mu^{2})=1/\mu^{2}$. The values $\xi=0.7$ and $\xi=1$ are
obtained from a linear extrapolation.

The form factors $\lambda_{i}(k,p)$, $i=1,2,3$, and $\tau_{i}(k,p)$,
$i=3,5,8$, are proportional to the ghost-dressing function $G(q^{2})$ defined
by the propagator,

$D^{ab}(q^{2})=-\,\delta^{ab}\,\frac{G(q^{2})}{q^{2}}\ ,$ (19)

and renormalized as $G(\mu^{2})=1$; $X_{0}^{\xi}(q^{2})$ is the leading form
factor of the quark-ghost scattering amplitude Marciano and Pagels (1978);
Davydychev et al. (2001), $H^{a\xi}(k,p)=H^{\xi}(k,p)t^{a}$, which can most
generally be decomposed as,

$\displaystyle H^{\xi}(k,p)$
$\displaystyle=\,X_{0}^{\xi}(k,p)\mathbbm{1}_{D}+iX_{1}^{\xi}(k,p)\,\gamma\cdot
k$ $\displaystyle+\,iX_{2}^{\xi}(k,p)\ \gamma\cdot
p+X_{3}^{\xi}(k,p)\,\left[\gamma\cdot k,\gamma\cdot p\,\right]\,.$ (20)

Calculated in one-loop dressed approximation Aguilar and Papavassiliou (2011);
Rojas et al. (2013), $X_{0}^{\xi}(k,p)$ can be projected out from the integral
equation for the quark-ghost scattering amplitude. We do so in the simplified
kinematic configuration $k=-p=q/2$, thus omitting an angular dependence
between $k$ and $p$, which is expressed by the following integral:

$\displaystyle X_{0}^{\xi}(q^{2})$
$\displaystyle=\,\tfrac{1}{4}\operatorname{Tr}_{\mathrm{CD}}\,H^{\xi}(q/2,-q/2)$
$\displaystyle=\,1+\frac{C_{A}}{4}\,g^{2}\\!\int^{\Lambda}\\!\\!\frac{d^{4}\ell}{(2\pi)^{4}}\
\Delta_{\mu\nu}^{\xi}(\ell)D(\ell+q)$
$\displaystyle\times\,\operatorname{Tr}_{D}\big{[}G_{\nu}\,S_{\xi}(\ell+\tfrac{q}{2})\,\Gamma_{\mu}^{\xi}(\ell+\tfrac{q}{2},-\tfrac{q}{2})\big{]}\
.$ (21)

Here, $C_{A}=3$ is the Casimir operator in the adjoint representation,

$G_{\nu}=i(\ell_{\nu}+q_{\nu})H_{1}-i\ell_{\nu}H_{2}\ ,$ (22)

is the dressed ghost-gluon vertex for which the form factor $H_{1}$ was
calculated in Ref. Dudal et al. (2012), $\ell$ is the gluon momentum exchanged
between the quark and the ghost and in this configuration $q$ coincides with
the gluon momentum in the DSE (II). The trace is over color and Dirac indices.
Beyond Landau gauge, the integral in Eq. (21) diverges since the gluon
propagator (2) is not transverse anymore and does therefore not project out
the $\ell_{\mu}$ terms of the ghost-gluon vertex (22). The divergence is
mildly logarithmic and given that we employ lattice-QCD propagators for
momenta $p\lesssim 8\,$GeV, no renormalization constant is necessary in the
numerical integration, even for large $\Lambda$. Nonetheless, we multiply the
left-hand side of Eq (21) with a renormalization factor $Z_{H}$ and impose
$X_{0}(\mu^{2})=1$.

We neglect the remaining $X_{i}^{\xi}(k,p)$ form factors as they are
suppressed with respect to $X_{0}^{\xi}(k,p)$ Aguilar et al. (2017a) in Landau
gauge. One may ask whether this approximation is justified in other gauges as
the form factors may vary in strength with the gauge parameter. This
dependence is depicted in Fig. 1 for $X_{0}^{\xi}(q^{2})$ and discussed in
more detail in Section III. It turns out that setting $X_{0}^{\xi}(q^{2})=1$
has only modest quantitative effects on $M_{\xi}(p^{2})$ and $Z_{\xi}(p^{2})$,
at least for $\xi\in[0,1]$. Therefore, we make bold to assume that our
observation is valid for the other form factors and will analyze their
behavior as a function of $\xi$ in a forthcoming study. We will see in Sec.
III that the contribution of $X_{0}(q^{2})$ to DCSB is small compared with the
effect of the transverse vertex.

Figure 1: Gauge-parameter dependence of the form factor $X_{0}^{\xi}(q^{2})$
associated with the quark-ghost scattering amplitude.  Figure 2: Mass function
$M_{\xi}(p^{2})$ and wave renormalization $Z_{\xi}(p^{2})$ as functions of the
gauge parameter $\xi$ obtained by solving the DSE (III.1) with the gluon and
ghost propagators in Eqs. (23) and (24), respectively, and the vertex
$\Gamma_{\mu}^{\xi}(k,p)$ of Eq. (6) with corresponding form factors
$\lambda_{i}^{\xi}(k,p)$ and $\tau_{i}^{\xi}(k,p)$.

The $Y_{i}^{(A,S)}$ functions represent the form-factor decomposition of the
Fourier transform of a four-point function in coordinate space in the TSTI.
The latter involves a line integral over a nonlocal vector vertex and a Wilson
line to preserve gauge invariance. As its matrix elements are rather
complicated expressions, we refer to the discussion in Ref. He (2009). We note
that the $Y_{i}^{(A,S)}$ functions are _a priori_ unknown and have been
constrained by us Albino et al. (2019) with the vertex ansatz in Ref. Bashir
et al. (2012b), i.e. we equated our transverse form factors (11)–(18) with
those of Ref. Bashir et al. (2012b) and derived expressions for the
$Y_{i}(k,p)$ that depend on $A(p^{2})$, $B(p^{2})$ and $X_{0}(p^{2})$. On the
other hand, if the $Y_{i}^{(A,S)}$ form factors were known exactly, so would
be the vertex. However, additional input, e.g. from the generalized LKFT that
describes the vertex itself, would be required for its full determination. We
remind that we do not solve the BSE of the quark-gluon vertex, hence we do not
truncate the vertex form factors. Therefore, the expressions in Eqs. (7) to
(18) represent a most complete form factor decomposition of the vertex that
satisfies SU(3) gauge symmetry, the $C$ and $P$ symmetries of the bare vertex
and which is not plagued by kinematic singularities in the limit $k\to p$.

In principle, one must also solve the gluon and ghost DSE which are coupled to
the gap equation (II). We deliberately abstain from doing so, as our interest
is in employing quenched gluon and ghost propagators in $R_{\xi}$ gauge from
lattice QCD Bicudo et al. (2015). For practical reasons in their numerical
implementation, we fit the gluon dressing functions with the parametrization
Dudal et al. (2018),

$\displaystyle\Delta_{\xi}(q^{2})=\frac{Z(q^{2}+M_{1}^{2})}{q^{4}+M_{2}^{2}q^{2}+M_{3}^{4}}\left[1+\omega\ln\left(\frac{q^{2}+M_{0}^{2}}{\Lambda_{\mathrm{QCD}}^{2}}\right)\right]^{-\gamma_{\mathrm{gl}}},$
(23)

where $\omega=11N_{c}\,\alpha_{s}(\mu)/12\pi$, $\Lambda_{\mathrm{QCD}}=0.425$
GeV and $\gamma_{\mathrm{gl}}=(13-3\xi)/22$ is the 1-loop anomalous gluon
dimension. The renormalization scale is $\mu=4.3$ GeV at which the strong
coupling is chosen to be $\alpha_{s}=0.29$ in Landau gauge Aguilar et al.
(2017b). This parametrization is motivated by the refined Gribov-Zwanziger
tree-level gluon propagator in the infrared domain and by the 1-loop
renormalization group behavior for large momenta, which amounts to a
renormalization-group improved Padé approximation. We collect the values $Z$,
$M_{0}$, $M_{1}$, $M_{2}$ and $M_{3}$ as a function of the gauge parameter
$\xi$ in Table 1. It turns out that all fit parameters depend linearly on
$\xi$ and we take this feature to our advantage to extrapolate the
parametrization to $\xi=1$.

Figure 3: Increasing and decreasing dependence of the mass function
$M_{\xi}(0)$ and wave renormalization $Z_{\xi}(0)$, respectively, on the gauge
parameter $\xi$.  Figure 4: The form factors $\lambda_{1}^{\xi}(-p,p)$ and
$\lambda_{3}^{\xi}(-p,p)$ as functions of the gauge parameter $\xi$, obtained
with Eqs. (7) and (9) in the symmetric limit $k=-p$. Note that the
dimensionless quantity $p\lambda_{3}^{\xi}$ is plotted.

The bare lattice data Duarte et al. (2016) for the ghost propagator is
parametrized with an analogous expression,

$\displaystyle G(q^{2})$
$\displaystyle=\frac{Z(q^{4}+M_{2}^{2}q^{2}+M_{1}^{4})}{q^{4}+M_{4}^{2}q^{2}+M_{3}^{4}}$
$\displaystyle\times\left[1+\omega\ln\left(\frac{q^{2}+\frac{m_{1}^{4}}{q^{2}+m_{0}^{2}}}{\Lambda_{\mathrm{QCD}}^{2}}\right)\right]^{\gamma_{\mathrm{gh}}},$
(24)

independent of $\xi$, as supported by a preliminary study in lattice QCD
Cucchieri et al. (2018). The anomalous ghost dimension is
$\gamma_{\mathrm{gh}}=-9/44$, while $\omega$, $\Lambda_{\mathrm{QCD}}$ and
$\mu$ are as in Eq. (23). A least-squares fit yields
$\chi^{2}/\mathrm{dof}=0.247$ and the parameters: $Z=5.068\pm 0.012$,
$M_{1}^{4}=19.281\pm 0.552\,\mathrm{GeV}^{4}$, $M_{2}^{2}=27.721\pm
0.696\,\mathrm{GeV}^{2}$, $M_{3}^{4}=7.695\pm 0.329\,\mathrm{GeV}^{4}$,
$M_{4}^{2}=24.340\pm 0.565\,\mathrm{GeV}^{2}$, $m^{2}_{0}=0.527\pm
1.263\,\mathrm{GeV}^{2}$, $m^{4}_{1}=0.018\pm 0.035\,\mathrm{GeV}^{4}$; this
fit must be normalized by a factor $N=1/4.706$ to ensure the renormalization
condition $G(4.3\,\mathrm{GeV})=1$.

##  III Gauge dependence of the quark propagator and invariance of the quark
condensate

With all the calculational tools and elements at hand in the gauge-parameter
interval between the Landau and Feynman gauges, we can evaluate the gauge
dependence of the quark’s mass and wave-renormalization functions as well as
that of the quark condensate.

###  III.1 Gauge dependent mass and wave-renormalization functions

Since the gauge propagators and the quark-gluon vertex are now determined we
can proceed to solve the DSE (II). After taking the color trace, and with the
renormalization procedure detailed in Ref. Albino et al. (2021), the DSE
becomes:

$\displaystyle S_{\xi}^{-1}(p)$ $\displaystyle=\,Z_{2}\,i\,\gamma\cdot
p+Z_{4}\,m(\mu)$ $\displaystyle+\,Z_{2}$
$\displaystyle\frac{16\pi\alpha_{s}^{\xi}}{3}\\!\int^{\Lambda}\\!\\!\frac{d^{4}k}{(2\pi)^{4}}\
\Delta^{\xi}_{\mu\nu}(q)\,\gamma_{\mu}\,S_{\xi}(k)\,\Gamma_{\nu}^{\xi}(k,p)\
.$ (25)

We plot the solutions of this DSE in Fig. 2 which clearly exposes that the
mass function $M_{\xi}(p^{2})$ increases with $\xi$ below the renormalization
point $\mu=4.3$ GeV, whereas the wave renormalization $Z_{\xi}(p^{2})$ is
suppressed. In fact, $M_{\xi}(0$) increases rather continuously with $\xi$
whereas $Z_{\xi}(0)$ slowly decreases as function of $\xi$, as can be read
from Fig. 3. We remark that the gauge dependence is threefold in Eq. (III.1),
namely $\xi$ enters directly via the longitudinal component and indirectly via
the dressing function of the transverse component, $\Delta_{\xi}(q^{2})$, of
the gluon propagator (2), but also via the gauge-dependent strong coupling
$\alpha_{s}^{\xi}$ for which we use the parametrization Aguilar et al.
(2017b):

$\alpha_{s}^{\xi}=0.29+0.098\xi-0.064\xi^{2}\ .$ (26)

The mass-function enhancement in Fig. 2 is due to the feedback of the gauge
dependent $A_{\xi}(p^{2})$ and $B_{\xi}(p^{2})$ functions in the dressed
quark-gluon vertex, namely Eqs. (7) to (18), as $M_{\xi}(p^{2})$ also
increases when we keep the coupling $\alpha_{s}=0.29$ constant for all values
of $\xi$. Moreover, $M_{\xi}(p^{2})$ slightly decreases with $\xi$ in the
rainbow truncation, $\Gamma^{\xi}_{\mu}(k,p)=\gamma_{\mu}$, while the DSE
solutions employing $\Gamma_{\mu}^{L\xi}(k,p)$ with Eqs. (7) to (10) show an
enhancement of the mass as a function of $\xi$. We emphasize that the
overwhelming contribution to mass generation is due to the massive terms in
the expressions for $\lambda_{i}^{\xi}(k,p)$ and $\tau_{i}^{\xi}(k,p)$. While
$X_{0}^{\xi}(q^{2})$ is gauge dependent, as becomes clear from Fig. 1 where an
enhancement of this form factor with increasing values of $\xi$ in the low-
momentum region is observed, setting $X_{0}^{\xi}(q^{2})=1,\xi\in[0,1]$, leads
to a 10% reduction of $M_{\xi=1}(0)$ in Feynman gauge and to almost no
variation of $M_{\xi=0}(p^{2})$ in Landau gauge. Neglecting
$X_{0}^{\xi}(q^{2})$ does not qualitatively alter our results, in particular
our conclusions about the quark condensate remain the same.

We also illustrate the gauge-variation impact on the strength of the quark-
gluon vertex in Fig. 4 with the form factors $\lambda_{1}^{\xi}(k,p)$ and
$\lambda_{3}^{\xi}(k,p)$ obtained in the symmetric limit $k=-p$ and with the
solutions for $M_{\xi}(p^{2})$ and $Z_{\xi}(p^{2})$ in Fig. 2. Both form
factors are enhanced, though this is more so the case for
$\lambda_{3}^{\xi}(-p,p)$. The latter is proportional to the mass function
$B_{\xi}(p^{2})$ whose variation with $\xi$ is more prominent.

###  III.2 Quark condensate and gauge invariance

As a practical application, we calculate the quark condensate which is an
order parameter for DCSB and was shown to be a manifestly gauge-invariant
quantity in any SU($N$) theory using the generalized LKFT (ABG)
transformations Aslam et al. (2016):

$-\langle\bar{q}q\rangle^{0}_{\xi}\,\equiv\,Z_{4}N_{c}\int^{\Lambda}\\!\frac{d^{4}k}{(2\pi)^{4}}\,\operatorname{tr}_{D}\left[S_{\xi}^{0}(k)\right]\,.$
(27)

To this end, we obtain $M_{\xi}(p^{2})$ and $Z_{\xi}(p^{2})$ in the limit
$m(\mu)\to 0$ with which we compute
$(-\langle\bar{q}q\rangle_{\xi}^{0})^{1/3}$ as a function of $\xi$. The result
is presented in Fig. 5, where the central value exhibits a moderate dependence
on $\xi$, i.e. a maximum deviation from the Landau-gauge value of 7% for
$\xi\gtrsim 0.3$ after which its dependence on $\xi$ appears to become milder.
The green-shaded band depicts an error estimate due to the statistical error
of $\pm 10$% of the gluon propagator. This is because the lattice calculations
of the gluon propagator in $R_{\xi}$ gauges were performed in smaller physical
volumes, $V\simeq(3.2\,\text{fm})^{4}$, and smaller gauge ensembles than for
the ghost propagator for which the fit reproduces the simulation results in a
large physical volume, $V\simeq(8.1\,\text{fm})^{4}$, and with an ensemble of
gauge configurations that is about three times larger. Therefore, the main
error source is due to the gluon propagator and we neglect a statistical error
associated with the ghost.

Our error estimate only enters via the transverse part of the gluon propagator
(2), yet with increasing values of the gauge parameter the contribution of the
longitudinal component is more important. This is the reason why the error
band is wider in Landau gauge and narrows towards Feynman gauge. We did not
include the uncertainty of the strong coupling (26) and the systematic error
due to our treatment of the nonlocal four-point functions in the TSTI.

Figure 5: Gauge dependence of the quark condensate. The horizontal pink-shaded
band indicates the admissible region of a gauge-independent chiral quark
condensate as implied by the LKFT in QCD

##  IV Final remarks

We have calculated the quark propagator’s dependence on the gauge parameter
$\xi$ in covariant gauges. It turns out that the gluon propagator
parameterization and $M_{\xi}(0)$ can be fitted with linear functions of $\xi$
within the gauge-parameter interval, $\xi\in[0,0.5]$, while $Z_{\xi}(0)$
decreases in a proportionate manner. Furthermore, the strength of the quark-
gluon interaction increases with $\xi$ which illustrates the gauge dependence
of the quark dynamics. Having said that, studying the local gauge
transformation of the nonperturbative quark propagator directly from its
solution in a given gauge is illuminating, as were similar QED3 studies Bashir
and Raya (2007) that explicitly demonstrated the gauge independence of the
fermion condensate. Both the Nielsen identities and the generalized LKFTs (ABG
and MDSDB) De Meerleer et al. (2020); Nielsen (1975); Aslam et al. (2016); De
Meerleer et al. (2018) formally provide the basis for this endeavor. However,
it is not a straightforward exercise in QCD, even at the perturbative level
Dall’Olio et al. (2021).

While Green functions are not physical observables, they are essential objects
for the elucidation and understanding of strong interactions. They enter each
and every hadron observable computed from QCD’s elementary degrees of freedom,
namely quarks and gluons (and ghosts within covariant gauges). How the
dependence on $\xi$ is washed out in the building of physical observables is a
nontrivial problem that constrains the truncated kernels of Bethe-Salpeter and
Faddeev equations and deserves further attention. The quark mass function
inserted in a bound-state equation, conspires with an intricate interaction
kernel involving quarks and gluons at all energy scales in a manner that
preserves axial vector Ward-Takahashi identity in addition to all other gauge
identities and transformations mentioned in the text in detail, yielding
gauge-independent physical observables. Any variation of these observables
with the gauge parameter is a measure of our departure from the full
implementation of the generalized LKFT and the Nielsen identities.

Our analysis of the DSE relies on the dressed vertex (6) described by form
factors $\lambda_{i}^{\xi}(k,p)$ and $\tau_{i}^{\xi}(k,p)$, which are plotted
for $\xi=0$ in Ref. El-Bennich et al. (2022) and show to be in reasonable
agreement with those obtained in lattice QCD simulations Kızılersü et al.
(2021). We are thus encouraged by our exploratory study on the gauge
dependence of the quark propagator, as they provide a valuable initial step in
the numerical representation of the generalized LKFT and Nielsen identities in
QCD despite the current limitations. However, further efforts must be made to
include all $X_{i}^{\xi}(k,p)$ form factors of the quark-ghost scattering
amplitude along with their angular dependence, to derive a true analytic
expression for the four-point function in the TSTI and to explore the
dependence on these hitherto ignored features in the current analysis.

## Acknowledgements

The work of A.B. is supported in part by the US Department of Energy (DOE)
Contract no. DE-AC05-06OR23177, under which Jefferson Science Associates, LLC
operates Jefferson Lab. It received funds from the Coordinación de la
Investigación Científica (CIC) of the University of Michoacán through grant
no. 4.10. B.E. receives financial support from the Brazilian agencies FAPESP,
grant no. 2018/20218-4, and CNPq, grant no. 428003/2018-4. F.E.S. is a CAPES-
PNPD postdoctoral fellow, grant no. 88882.314890/2013-01. J.R.S. was supported
by a CAPES PhD fellowship. O.O. acknowledges support from FCT, Portugal,
project nos. UID/FIS/04564/2019 and UID/FIS/04564/2020. This work is part of
the network project INCT-Física Nuclear e Aplicações, no. 464898/2014-5.

## References

  * Ward (1950) J. C. Ward, Phys. Rev. 78, 182 (1950).
  * Takahashi (1957) Y. Takahashi, Nuovo Cim. 6, 371 (1957).
  * Slavnov (1972) A. A. Slavnov, Theor. Math. Phys. 10, 99 (1972).
  * Taylor (1971) J. C. Taylor, Nucl. Phys. B 33, 436 (1971).
  * Ball and Chiu (1980) J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2542 (1980).
  * Aguilar et al. (2017a) A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 96, 014029 (2017a), eprint 1610.06158.
  * Bermudez et al. (2017) R. Bermudez, L. Albino, L. X. Gutiérrez-Guerrero, M. E. Tejeda-Yeomans, and A. Bashir, Phys. Rev. D 95, 034041 (2017), eprint 1702.04437.
  * Takahashi (1985) Y. Takahashi, in _Positano Symp.1985:0019_ (1985), p. 0019.
  * Kondo (1997) K.-I. Kondo, Int. J. Mod. Phys. A12, 5651 (1997), eprint hep-th/9608100.
  * He et al. (2000) H.-X. He, F. C. Khanna, and Y. Takahashi, Phys. Lett. B480, 222 (2000).
  * He (2006) H.-X. He, Commun. Theor. Phys. 46, 109 (2006).
  * He (2007) H.-X. He, Int. J. Mod. Phys. A22, 2119 (2007).
  * He (2009) H.-x. He, Phys. Rev. D 80, 016004 (2009), eprint 0906.2834.
  * Qin et al. (2013) S.-X. Qin, L. Chang, Y.-X. Liu, C. D. Roberts, and S. M. Schmidt, Phys. Lett. B 722, 384 (2013), eprint 1302.3276.
  * Qin et al. (2014) S.-X. Qin, C. D. Roberts, and S. M. Schmidt, Phys. Lett. B733, 202 (2014), eprint 1402.1176.
  * Albino et al. (2019) L. Albino, A. Bashir, L. X. G. Guerrero, B. E. Bennich, and E. Rojas, Phys. Rev. D 100, 054028 (2019), eprint 1812.02280.
  * Albino et al. (2021) L. Albino, A. Bashir, B. El-Bennich, E. Rojas, F. E. Serna, and R. C. da Silveira, JHEP 11, 196 (2021), eprint 2108.06204.
  * Landau and Khalatnikov (1955) L. D. Landau and I. M. Khalatnikov, Zh. Eksp. Teor. Fiz. 29, 89 (1955).
  * Fradkin (1955) E. S. Fradkin, Zh. Eksp. Teor. Fiz. 29, 258 (1955).
  * Aslam et al. (2016) M. J. Aslam, A. Bashir, and L. X. Gutierrez-Guerrero, Phys. Rev. D 93, 076001 (2016), eprint 1505.02645.
  * De Meerleer et al. (2018) T. De Meerleer, D. Dudal, S. P. Sorella, P. Dall’Olio, and A. Bashir, Phys. Rev. D 97, 074017 (2018), eprint 1801.01703.
  * De Meerleer et al. (2020) T. De Meerleer, D. Dudal, S. P. Sorella, P. Dall’Olio, and A. Bashir, Phys. Rev. D 101, 085005 (2020).
  * Nielsen (1975) N. Nielsen, Nuclear Physics B 101, 173 (1975).
  * Bashir et al. (2011) A. Bashir, A. Raya, and S. Sanchez-Madrigal, Phys. Rev. D 84, 036013 (2011), eprint 1108.4748.
  * Bashir et al. (1999) A. Bashir, A. Kizilersu, and M. R. Pennington (1999), eprint hep-ph/9907418.
  * Bashir et al. (2000) A. Bashir, A. Kizilersu, and M. R. Pennington, Phys. Rev. D 62, 085002 (2000), eprint hep-th/0010210.
  * Bashir and Raya (2005) A. Bashir and A. Raya, Nucl. Phys. B 709, 307 (2005), eprint hep-ph/0405142.
  * Bashir and Raya (2007) A. Bashir and A. Raya, Few Body Syst. 41, 185 (2007), eprint hep-ph/0511291.
  * Bicudo et al. (2015) P. Bicudo, D. Binosi, N. Cardoso, O. Oliveira, and P. J. Silva, Phys. Rev. D 92, 114514 (2015), eprint 1505.05897.
  * Bashir et al. (2012a) A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y.-X. Liu, C. D. Roberts, and P. C. Tandy, Commun. Theor. Phys. 58, 79 (2012a), eprint 1201.3366.
  * Huber (2015) M. Q. Huber, Phys. Rev. D 91, 085018 (2015), eprint 1502.04057.
  * Napetschnig et al. (2021) M. Napetschnig, R. Alkofer, M. Q. Huber, and J. M. Pawlowski, Phys. Rev. D 104, 054003 (2021), eprint 2106.12559.
  * Bowman et al. (2004) P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, and A. G. Williams, Phys. Rev. D 70, 034509 (2004), eprint hep-lat/0402032.
  * Oliveira et al. (2019) O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Phys. Rev. D 99, 094506 (2019), eprint 1809.02541.
  * Rojas et al. (2013) E. Rojas, J. P. B. C. de Melo, B. El-Bennich, O. Oliveira, and T. Frederico, JHEP 10, 193 (2013), eprint 1306.3022.
  * Serna et al. (2019) F. E. Serna, C. Chen, and B. El-Bennich, Phys. Rev. D 99, 094027 (2019), eprint 1812.01096.
  * El-Bennich et al. (2022) B. El-Bennich, F. E. Serna, R. C. da Silveira, L. A. F. Rangel, A. Bashir, and E. Rojas, Rev. Mex. Fis. Suppl. 3, 0308092 (2022), eprint 2201.04144.
  * Marciano and Pagels (1978) W. J. Marciano and H. Pagels, Phys. Rept. 36, 137 (1978).
  * Davydychev et al. (2001) A. I. Davydychev, P. Osland, and L. Saks, Phys. Rev. D 63, 014022 (2001), eprint hep-ph/0008171.
  * Aguilar and Papavassiliou (2011) A. C. Aguilar and J. Papavassiliou, Phys. Rev. D 83, 014013 (2011), eprint 1010.5815.
  * Dudal et al. (2012) D. Dudal, O. Oliveira, and J. Rodriguez-Quintero, Phys. Rev. D 86, 105005 (2012), eprint 1207.5118.
  * Bashir et al. (2012b) A. Bashir, R. Bermudez, L. Chang, and C. D. Roberts, Phys. Rev. C 85, 045205 (2012b), eprint 1112.4847.
  * Dudal et al. (2018) D. Dudal, O. Oliveira, and P. J. Silva, Annals Phys. 397, 351 (2018), eprint 1803.02281.
  * Aguilar et al. (2017b) A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D 95, 034017 (2017b), eprint 1611.02096.
  * Duarte et al. (2016) A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 014502 (2016), eprint 1605.00594.
  * Cucchieri et al. (2018) A. Cucchieri, D. Dudal, T. Mendes, O. Oliveira, M. Roelfs, and P. J. Silva, PoS LATTICE2018, 252 (2018), eprint 1811.11521.
  * Dall’Olio et al. (2021) P. Dall’Olio, T. De Meerleer, D. Dudal, S. P. Sorella, and A. Bashir, Nucl. Phys. B 973, 115606 (2021), eprint 2107.03910.
  * Kızılersü et al. (2021) A. Kızılersü, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Phys. Rev. D 103, 114515 (2021), eprint 2103.02945.