# The analytic structure of the Landau gauge quark propagator from Padé
analysis

Alexandre F. Falcão${}^{1,2\,}$<EMAIL_ADDRESS>Orlando
Oliveira${}^{2}\,$<EMAIL_ADDRESS>
1 Department of Physics and Technology, University of Bergen, 5007 Bergen,
Norway

2CFisUC, Departament of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal

###### Abstract

The analytic structure of the 2 flavour full QCD lattice Landau gauge quark
propagator is investigated with Padé approximants applied to its vector and
scalar form factors. No poles at complex momentum are observed for the
propagator. Moreover, there is clear evidence of a pole at real on-axis
negative Euclidean momentum, i.e. for Minkowski type of momentum. This pole
occurs at Euclidean momenta $p^{2}\sim-300$ MeV and it reproduces typical
quark mass values used in phenomenological effective quark models. The Padé
approximant analysis also gives hints on the presence of a branch cut. Our
results also show a clear correlation between the position of this pole,
understood as an effective quark mass, and the pion mass that is compatible
with PCAC. Slightly differences between the poles for the two quark form
factors are observed which can be viewed either as a limitation of the method
or as a suggestion that the quark propagator has no spectral representation.

###### Contents

  1. I Introduction and Motivation
  2. II Padé approximants
  3. III Resumé of the quenched analysis
    1. III.1 The Ghost Propagator
    2. III.2 The Gluon Propagator
  4. IV The Quark Propagator
  5. V Summary and Conclusions

##  I Introduction and Motivation

The theory that describes the color interaction between quarks and gluons, i.e
Quantum Chromodynamics (QCD) (see Alkofer:2000wg ; Fischer:2006ub ;
Binosi:2009qm  for reviews), should also explain hadron phenomenology.
However, so far, in what concerns the experimental evidences for free quarks
and gluons Perl:2009zz ; Tanabashi:2018oca  there are no positive results.
This lack of experimental evidences suggests that quarks and gluons can only
appear as components of the observed particles, the hadrons, that are
understood as color singlet objects. This negative results are translated into
the confinement hypothesis that still requires to be prooved. Despite the
success of QCD as a fundamental theory, the bridge linking quarks and gluons
to the observed hadrons still needs to be paved. From the theoretical side,
the understanding of quarks and gluons as components of hadrons requires the
knowledge of the fundamental QCD Green functions such as the propagators and
vertices, as e.g. the quark-gluon vertex and the three-gluon vertex.

In the functional continuum approach to QCD, the solution of the Dyson-
Schwinger equations for the quark, the ghost and the gluon propagators are
used to feed Bethe-Salpeter equations, Faddeev equations, etc. to access the
properties of the physical hadronic states. The computation of the solutions
of these equations requires the knowledge of the propagators in the full
complex momentum plane. Therefore, the analytic structure of the propagators
have importante roles. Indeed, being able to identify the poles and branch
cuts of the propagators is crucial to the bound state problem but also for the
comprehension of the confinement mechanism, to the dynamical properties of
quarks and gluons and, ultimately to the understanding of QCD itself.

For the gluon and ghost sectors of QCD, the past twenty years saw an effort to
have reliable descriptions and interpretations for the propagators. A recent
review on the gluon propagator can be found in Papavassiliou:2022wrb . See,
also, Leinweber:1998uu ; Becirevic:1999uc ; Becirevic:1999hj ;
Cucchieri:2007md ; Bogolubsky:2009dc ; Dudal:2010tf ; Cucchieri:2011ig ;
Maas:2011ez ; Dudal:2018cli ; Li:2019hyv ; Catumba:2021hcx  and references
therein for lattice calculations and vonSmekal:1997ohs ; Alkofer:2003jj ;
Boucaud:2007va ; Boucaud:2007hy ; Huber:2007kc ; Dudal:2008rm ; Aguilar:2008xm
; Rodriguez-Quintero:2010qad ; Aguilar:2010zx ; Strauss:2012zz ;
Eichmann:2021zuv ; Pelaez:2014mxa ; Aguilar:2019uob ; Reinosa:2020skx ;
Fischer:2020xnb ; Pelaez:2021tpq ; Dudal:2022nnu  and references therein for
continuum approaches. On the other hand, the ghost propagator has also been
thoroughly studied, as can be found in e.g. vonSmekal:1997ohs ; Alkofer:2003jr
; Aguilar:2008xm ; Bogolubsky:2009dc ; Duarte:2016iko ; Boucaud:2017ksi ;
Duarte:2017wte  and references therein. Similarly, one can find in the
literature several studies for the quark propagator, see e.g Krein:1990sf ;
Bowman:2002bm ; Boucaud:2003dx ; Fischer:2005nf ; Furui:2006rx ; Furui:2006ks
; Ayala:2012pb ; Oliveira:2018lln ; Solis:2019fzm ; Lessa:2022wqc , that
include continuum and lattice QCD calculations. From the good agreement
obtained by several groups that use different techniques, it seems reasonable
to claim a fair understanding of the QCD two point correlation functions.
However, most non-perturbative computations are done using the Euclidean
formulation of the theory, which limits the determination of the analytical
structure of the propagators.

The quark propagator has been investigated by several groups Maris:1994ux ;
Maris:1995ns ; Burden:1997ja ; Bhagwat:2003vw ; Alkofer:2003jj ;
Frederico:2019noo ; Hayashi:2021nnj ; Hayashi:2021jju , with the Dyson-
Schwinger studies suggesting the presence of poles at complex momenta.
However, it is difficult to disentangle if the complex poles and analytic
structure result directly from the truncations and/or the parameterization of
the vertices used to solve the Dyson-Schwinger equations themselves.
Therefore, it is important to investigate the analytical structure of the QCD
propagators with different methods, different sets of data, and to confront
the outcome of different techniques.

The main goal of the current work is to address the analytic structure of the
quark propagator by studying a subset of the Landau gauge lattice data
generated in Oliveira:2018lln , by applying the techniques based on Padé
approximants considered in Falcao:2020vyr  for the pure Yang-Mills gluon and
ghost propagators.

In this previous work Falcao:2020vyr , the study of the analytic structure of
the pure Yang-Mills theory gluon and ghost propagators used only the absolute
value of the residua and their relative importance. Herein, we take the
opportunity to complement this information and look also at the sign of the
relevant residua. The sign of the residua provides a naive and clear
interpretation of the corresponding poles. In this respect, for real on-axis
momenta, see Fig. 1, some of the ghost residua are negative and the poles
cannot represent physical particles. For the gluon propagator and for real on-
axis momenta, the interpretation of the residua is not so clear, see Fig. 2,
and no conclusions can be established firmly. As discussed in Falcao:2020vyr ,
the Padé analysis for the gluon propagator identifies clearly a pair of
complex conjugate poles that are not seen in the ghost propagator.

The analysis of quark propagator data with Padé approximants discussed below
shows an analytical structure that differs from the analytic structure of the
gluon and the ghost propagators. Contrary to the Dyson-Schwinger studies for
complex momenta, no evidence is found for poles at complex momenta. On the
other hand, the Padé analysis gives indications of a possible branch cut for
real on-axis momenta. Furthermore, the Padé approximant description of the
lattice data finds poles at real on-axis negative Euclidean momenta around
$p^{2}\sim-300$ MeV, i.e. for Minkowski momenta, that are correlated with the
pion mass. Indeed, it is observed that the correlation between the quark mass,
understood as the pole mass at the Minkowski momenta, and the pion mass is
compatible with the predictions of PCAC. Moreover, the residua of these poles
for time-like Euclidean momenta is positive. This seems to indicate that the
quark propagator has a particle-like pole at Minkowski momenta, at momenta
that approximately reproduce typical values for effective quark masses, and
further structures, that we are not able to disentangle, that must prevent the
quark from being a physical particle. Recall that, for example, the Dyson-
Schwinger studies for the quark propagator observe positivity violation.
Positivity violation implies that if a Kállën-Lehmann representation of the
propagator exists, then its spectral representation is not positive defined. A
spectral function that is not positive defined is not compatible with a
quantum mechanical probabilistic interpretation and single quark state cannot
belong to the set of physical states of QCD. Naively one could think that the
presence of the pole at Minkowsky type of momenta with a positive residua
implies that such type of pole would correspond to a physical particle but
this it not necessarily true. A counter example can be find in Solis:2019fzm ,
where the authors provide a quark propagator built from a toy model that have
similar properties as those just described but do not represent a physical
particle. The pole at the Minkowsky momenta occurs for momenta that are within
typicall values for the quark mass used in effective quark models. In this
sense, the results found herein seem to give some support to the effective
quark models and to the phenomenological values that they use for the quark
mass. A clarification of these questions calls for an understanding of what
are the signs of confinement in the quark propagator, a complex problem that
clear goes beyond the type of answers that we are able to provide in the
current work.

The manuscript is organised as follows. In Sec. II the Padé approximant theory
to address the analytical structure of the propagators is reviewed. Then, the
results for the pure Yang-Mills theory for the ghost and the gluon propagators
are revisited in Sec. III.1 and Sec. III.2, respectively. The analysis of
quark propagator data with Padé approximants is performed in Sec. IV. Finally,
in Sec. V we summarise the results and conclude.

##  II Padé approximants

The Padé approximant to a function $f(x)$ defined in the interval $x\in[a,b]$
looks for the best approximation to $f(x)$ written as a ratio of polynomials
in $x$. The polynomials associate zeros and poles to $f(x)$, corresponding to
the zeros of the numerator and denominator, respectively. The location of
these zeros and poles depends on the degree of the polynomials considered and
changing the degrees of the numerator and denominator polynomals also changes
the position of the zeros and the poles. However, in general, there is a
subset of zeros and poles whose position does not depend on the Padé
approximant used. These zeros and poles, which are independent of the Padé
approximant, can be associated with the analytic structure of $f(x)$. The
remaining zeros and poles are artefacts of the approximation. Moreover, a
stable structure of in the complex plane formed by a sequence of alternating
poles and zeros of the approximant can represent a branch cut.

For particle propagators in Quantum Field Theory represented, generally, as
$D(p^{2})$, the Padé approximant of order $[\,M\,|\,N\,]$ to the propagator is
given by

$D(p^{2})\approx P^{M}_{N}(p^{2})=\frac{Q_{M}(p^{2})}{R_{N}(p^{2})}\ ,$ (1)

with

$\displaystyle
Q_{M}(p^{2})=q_{0}+\cdots+q_{M}\,\left(p^{2}\right)^{M}\qquad\mbox{ and
}\qquad R_{N}(p^{2})=1+\cdots+r_{N}\,\left(p^{2}\right)^{N}$ (2)

and where, by convention, the coefficient of the lowest order term in the
denominator is set to one. The problem of building a Padé approximant to
$f(x)$ is reduced to find the set of coefficients
$q_{0},\,q_{1},\dots,\,r_{1},\dots$ that define the best approximation to
$D(p^{2})$ according to Eq. (1). Pommerenke’s Theorem Pomm73  states that for
a meromorphic function $f(z)$, the Padé sequences $[\,M\,|\,M+k\,]$, with
fixed $k$, converge to $f(z)$ in any compact set of the complex plane. This
theorem gives support to use the Padé approximant sequences in the
investigation of the analytic structure of particle propagators.

In the sequences of Padé approximants $[\,M\,|\,M+k\,]$, single poles of
$f(z)$ appear as stable poles for sufficiently large values of $M$. On the
other hand, the Froissart doublets Baker75 ; BaMo96 ; BeOr99 ;
SanzCillero:2010mp ; Queralt:2010sv  are those poles that depend strongly on
the degree of the polynomials or have nearby zeros. The presence of the nearby
zeros for the Froissart doublets result in residua whose absolute value is
small. The absolute value of the residua can then be used to distinguish the
meaningful poles from those that are artefacts of the method. The doublets
that appear at sufficiently large values of $M$ are artefacts associated with
the use of ratio of polynomials to represent $f(x)$.

It is common practice to use diagonal $[\,M\,|\,M\,]$ or near diagonal
$[\,M\,|\,M\pm 1\,]$ Padé approximants in the representation of $f(x)$. In our
analysis of the quark propagator, motivated by the studied of the pure Yang-
Mills gluon and ghost propagators Falcao:2020vyr , and also to be closer to
the perturbative description of the propagator, that should be recovered at
high momenta, only the sequence $[\,N\,|\,N+1\,]$ will be considered.

For the determination of the coefficients of the polynomials that define a
given Padé approximant we look at the absolute minimum of the objective
function

$\chi^{2}=\sum_{j=1}^{N}\left(\frac{D(p^{2}_{j})-D_{Lat}(p^{2}_{j})}{\sigma(p^{2}_{j})}\right)^{2}\
,$ (3)

where the sum is over the data points obtained with lattice QCD simulations,
$D(p^{2})=P^{N}_{N+1}(p^{2})$, $D_{Lat}(p^{2})$ are the data points for the
given function and $\sigma(p^{2})$ are the associated statistical error with
$D_{Lat}(p^{2})$. In the following we do not take into consideration
correlations between the momenta. In this way, the numerical problem becomes a
non-linear global optimisation problem that is solved with the help of the
global optimisation methods available within Mathematica Math  software
package, i.e. the Mathematica implementations of the differential evolution
(DE) method and of the simulated annealing (SA) method. We recall the reader
that global optimization problems are non-trivial per se and the comparison of
the two methods is important to have confidence in the results. Furthermore,
we report that in Falcao:2020vyr  the authors performed a series of tests to
understand the numerical performance of the numerical method before studying
the analytic structure of the gluon and ghost pure Yang-Mills propagators. We
refer the interested reader to this work for further details concerning
numerical experiments, performance and limitations of the method.

Figure 1: Padé estimations for the ghost propagator residua for real on-axis
momenta as a function of $p^{2}$ (left) and a function of the degree of the
Padé approximant (right).

##  III Resumé of the quenched analysis

As a warm up to the analysis of the Landau gauge full QCD quark propagator, we
revisit the results of Padé approximants analysis for the analytic structure
of the pure gauge gluon and ghost propagators performed in Falcao:2020vyr ,
taking the opportunity to comment on the value of the residua.

Figure 2: Padé estimations for the gluon propagator residua for real on-axis
momenta for the simulations using the two largest lattice volumes considered
in Falcao:2020vyr .

###  III.1 The Ghost Propagator

The analysis of the Landau gauge ghost propagator analytic structure concluded
for the absence of poles at complex momenta, for the presence of a pole at
zero momenta and found possible evidences for a branch cut at time-like
(Minkowski) momenta. The clear identification of the branch points is not well
established and it requires further studies, that call for simulations with
larger ensembles of gauge configurations and/or larger physical volumes. The
analysis of the analytic structure of the propagator in Falcao:2020vyr  was
centred on the absolute value and on the relative importance of the absolute
value of the residua.

The Padé estimations of the residua $Z$ for the ghost propagator for real on-
axis monenta are reported in Fig. 1. For real on-axis momenta, the residua are
all real numbers as required if the propagator is a real function of $p^{2}$.
The data shows clearly that the largest residua has $Z\sim 3$ and is
associated with the pole at zero momentum. This pole appears for all the $N$
considered in the analysis. The second highest residuum in absolute value is
negative and, therefore, cannot represent a physical particle. Furthermore, it
is associated with a negative $p^{2}$ that corresponds to Minkowsky or time-
like momenta. As seen in Fig. 1, there is a pole at $p^{2}\sim-0.5$ GeV2 that
both methods identify only for small $N$, and another pole at $p^{2}=-1$ GeV2
also with a negative residuum that is identified for $N\geqslant 11$ and whose
position seems to be independent of the degree of the approximant111 For
$N=11$ the differential evolution method identifies a pole at $p^{2}=-0.992$
GeV2 with a residuum of $Z=-0.4360$, while the simulated annenaling sees a
pole at $p^{2}=-0.984$ GeV2 with $Z=-0.4307$.. As discussed in Falcao:2020vyr
, see their Fig. 18, the pole associated with $p^{2}\sim-0.5$ and $-1$ Gev2
have an associated close by zero of the numerator that suggests the presence
of a branch point for Minkowski momenta. As observed in the above cited work,
there is another pole at smaller momentum, that is stable under variation of
degree of the Padé approximant, whose residua is negative and with no close by
zero of the polynomial approximation222For the differential evolution this
pole appears firstly for $N=6$ and for $p^{2}=-0.096$ GeV2 with a $Z=-1.006$,
while the simulated annealing returns for $N=6$ a $p^{2}=-0.042$ Gev2 for a
$Z=-0.548$.. Its residuum being real and negative, the pole cannot represent a
physical particle. The remaining observed poles have much smaller $|Z|$, their
position is not stable against the variation of $N$ and, therefore, we do not
consider them as having a physical meaning. We call the reader’s attention
that the analysis of the infrared region ($p^{2}<1$ GeV2) for the same ghost
propagator data performed in Li:2021wol , although in a particular theoretical
framework, also encountered a dominant pole at $p^{2}=0$ with a positive
residuum and a second pole at $p^{2}\sim-0.3$ GeV2 with a negative residuum
but with smaller absolute value.

###  III.2 The Gluon Propagator

The Padé analysis of the Landau gauge lattice gluon propagator identifies a
pair of complex conjugate poles at $p^{2}\sim-0.28\pm\,i\,0.43$ GeV2. No poles
or branch cuts are clearly identified. Indeed, as discussed in Falcao:2020vyr
, there are indications of a possible branch cut but, as for the ghost
propagator, no firm conclusions were drawn.

The residua for the gluon propagator estimated using Padé approximants for on-
axis momenta are all real numbers and are reported in Fig. 2 for the two
larger lattice volumes considered previously and using the two global
optimization methods, i.e. differential evolution and simulated annealing. As
can be observed in Fig. 2 there seems to be stable poles at $p^{2}\sim\pm 1$
GeV2 but their residua are small in absolute value and, for some of the $N$,
they have close by zeros; see Fig. 14 in Falcao:2020vyr . These two features
prevent us to assign a meaning to these poles besides being possible branch
points. Note also that, for on-axis momenta, the two global optimization
methods identify structures in the same range of momenta but their results are
not fully consistent.

##  IV The Quark Propagator

Let us now look at the analysis of the Landau gauge full quark propagator for
a subset of the data published in Oliveira:2018lln . In the following, we will
investigate the results of the ensembles simulated at $\beta=5.29$, that
corresponds to a lattice spacing of $a=0.071$ fm, and from the available data
only the data for pion masses of 422 MeV, 290 MeV and 150 MeV. The first
simulation was performed on an asymmetric $32^{3}\times 64$ lattice, while the
last two simulations used a $64^{4}$ lattice. The simulation with an almost
physical pion mass uses about half of the gauge configurations of the other
two cases and, therefore, the associated statistical errors are larger.

Figure 3: Lattice Landau gauge quark propagator data for the vector (left) and
scalar (right) form factors. See text for definitions.

Figure 4: Minimum of the $\chi^{2}/d.o.f.$ achieved by the global optimization
methods for each of the lattice data sets. From top to bottom, the first line
of plots refers to the $M_{\pi}=150$ MeV, the second line to $M_{\pi}=290$ MeV
and the bottom line to $M_{\pi}=422$ MeV. In each line, the left plot is the
outcome of the minimization of the vector form factor, while the right plot
refers to the minima for the scalar form factor.

Figure 5: Poles, zeros and absolute value of the residua for complex valued
momenta for the simulation with a $M_{\pi}=290$ MeV. On the upper plot the
data refers to the vector form factor, while the lower plot refers to the
analysis of the scalar form factor. Here the poles with $|Z|<0.5$ are not
shown.

Figure 6: Poles, zeros and residua for real on-axis momenta for $M_{\pi}=422$
MeV. The upper plot reports the results from the analysis of the vector form
factor, while the lower plot reports the results of the analysis of the scalar
part form factor.

The analysis of the Landau gauge lattice quark propagator assumes that this
two point correlation function is color diagonal and that its Lorentz-Dirac
structure reads, in momentum space,

$S(p)=Z(p^{2})\frac{\not{p}+M(p^{2})}{p^{2}+M^{2}(p^{2})}=\frac{Z(p^{2})}{D(p^{2})}\,\Big{(}\not{p}+M(p^{2})\Big{)}$
(4)

where $D(p^{2})=p^{2}+M^{2}(p^{2})$. In the following, we name

$\frac{Z(p^{2})}{D(p^{2})}\qquad\mbox{ and
}\qquad\frac{Z(p^{2})\,M^{2}(p^{2})}{D(p^{2})}$ (5)

as vector and scalar form factors, respectively. The lattice data for the
propagator for the various ensembles can be seen in Fig, 3. The vector and
scalar form factors were rebuilt from the original lattice data assuming
Gaussian distributions for the propagation of errors. As see in Fig. 3, the
vector and scalar form factors are enhanced at low momenta and the enhancement
increases as the quark mass, or the pion mass, decreases. However, for the
smaller quark mass, the scalar form factor drops faster as the momentum
increases, with the form factor becoming smaller than the two heavier pion
masses considered for momenta $p\gtrsim 0.5$ GeV. The comparison of the left
and right plots in Fig. 3 reveals a complex situation that is the result of
the dependence of the quark wave function $Z(p^{2})$ and running quark mass
$M(p^{2})$ with the quark mass, i.e. with the pion mass.

The lattice data is to be described by a Padé approximant where the
coefficients of the polynomials are determined by minimizing the corresponding
$\tilde{\chi}^{2}=\chi^{2}/d.o.f.$ The quality of the minimization can be
measured by the value of $\tilde{\chi}^{2}=\chi^{2}/d.o.f.$ at the minimum.
These values, for the various sets of data, for the two minimizating
algorithms and for the various $N$ are reported in Fig. 4. In general, the
values of the $\tilde{\chi}^{2}$ obtained for the Padé approximants are able
to reproduce well the quark propagator form factors, with the
$\tilde{\chi}^{2}$ at the minimum associated with the simulated annealing
method performing slightly worst than for the differential evolution method.
The exception is the scalar data for the heaviest pion mass $M_{\pi}=422$ MeV,
whose $\tilde{\chi}^{2}$ takes values well above the acceptable. In the
following, we will disregard the data coming from the analysis of the scalar
form factor associated with $M_{\pi}=422$ MeV.

The Padé study of the poles, the zeros and the residua for complex momenta of
the Padé approximants shows no stable structures for momenta $|p^{2}|\leqslant
2$ Gev2. The same was observed for the pure gauge ghost propagator in
Falcao:2020vyr . This is illustrated in Fig. 5 for the simulation performed
with $M_{\pi}=290$ MeV. Similar plots can be show for the other two data sets.
We take these results as an indication that the quark propagator has no poles
for complex momenta. This configures a quite different nature for the analytic
structure of the fundamental QCD propagators. If the pure gauge gluon
propagator has poles at complex momenta as suggest in Falcao:2020vyr ;
Binosi:2019ecz  and also by the good agreement found between the predictions
based on the use of the Gribov-Zwanziger actions Zwanziger:1989mf ;
Dudal:2008sp ; Dudal:2007cw  and the lattice data for the infrared region
Dudal:2010tf ; Oliveira:2012eh ; Cucchieri:2011ig ; Cucchieri:2016jwg ;
Dudal:2018cli , the Padé analysis suggests that the pure gauge ghost
propagator and the quark propagator are void of poles for complex momenta. If
this is the case, it remains to be understood the mechanism that in QCD
distinguishes the bosonic degrees of freedom from those whose nature is of the
Grassmann type.

Figure 7: The same as Fig. 6 but for the data associated with $M_{\pi}=290$
MeV.

Figure 8: The same as Figs 6 and 7 but for the data associated with
$M_{\pi}=150$ MeV.

Let us turn now our attention to on-axis momenta. The outcome of the analysis
of the form factors for $M_{\pi}=422$ MeV can be seen in Fig. 6. In what
concerns the vector form factor, the analysis of the results obtained with the
differential evolution method does not have a clear interpretation. On the
other hand, the results coming from the simulated annealing method reveal a
pole at time-like momenta $p^{2}\sim-0.35$ GeV2 whose position is essentially
independent of the degree of Padé approximant considered. The results for the
smaller $N$ where this pole appears, give a pole position at $p^{2}=-0.275$
GeV2 with a residuum $Z=1.535$ when $N=4$ and a pole at $p^{2}=-0.227$ GeV2
with $Z=1.266$ for $N=5$. The corresponding results obtained with the
differential evolution method are a pole $p^{2}=-0.164$ GeV2 and a residuum of
$Z=1.087$ for $N=4$ and pole at $p^{2}=-0.293$ GeV2 with $Z=2.325$ for $N=5$
are in the same ballpark as those obtained with the simulated annealing
method. From the data it is difficult to identify any other structure as, for
example, a possible branch point. As already mentioned, the results coming
from the analysis of the scalar form factor are not so faithful. However,
looking at the bottom plot of Fig. 6, it qualitatively reproduces the results
of the vector form factor with the poles at slightly smaller values of
$|p^{2}|$. This may suggest that the poles of the quark propagator do not have
to occur at the same momentum for its vectorial and scalar form factors. This
contrasts with the intuitive picture that is built from the analysis of the
free fermionic propagators as, for example, is found in QED or other theories
where fermions appear as free particles.

The results of the analysis of the lattice data for the simulation with an
$M_{\pi}=290$ MeV are reported in Fig. 7. For the vector form factor (upper
plot), the outcome of the simulated annealing method is again easier to
understand. It shows a pole at time-like momenta, that is also observed with
the differential evolution method. For the smallest values of $N$, the
simulated annealing returns a pole at $p^{2}=-0.139$ GeV2 and a residuum of
$Z=1.174$ for $N=2$ and $p^{2}=-0.148$ GeV2 with $Z=1.201$ for $N=3$, while
the results of the differential evolution method return a pole at
$p^{2}=-0.139$ GeV2 with $Z=1.174$ for $N=2$ and $p^{2}=-0.127$ GeV2 with
$Z=1.203$ for $N=4$. No further clear poles or possible branch points are
identified from the analysis of the vector form factor. In what concerns the
results for the scalar form factors, the two optimization methods give similar
results that sugest a branch point at $p^{2}\sim-1$ GeV2 where there is a
proliferation of poles that have close zeros of the Padé approximants. A
precise determination of a possible branch point would demand for less
scattered data that, in principle, is possible to achieve with higher
statistical simulations. The data for the scalar form factor shows a clear
pole whose location is at $p^{2}=-0.119$ GeV2 with a residuum $Z=0.299$ for
$N=8$ and $p^{2}=-0.223$ GeV2 with $Z=0.740$ for $N=10$ according to the
differential evolution method, and at $p^{2}=-0.205$ GeV2 with a residuum
$Z=0.727$ for $N=4$ and $p^{2}=-0.182$ GeV2 with $Z=0.623$ for $N=6$ by the
simulated annealing method. These results suggest that the pole of the vector
form factor is at $p^{2}=-0.2$ GeV2 and its residuum $Z=1.2$, while the pole
of the scalar form factor occurs at $p^{2}=-0.2$ GeV and its residuum is
$Z=0.65$. Recall that as for the largest pion mass, the pole residuum is
positive defined and the pole of the vector and scalar form factors occur for
the same time-like momentum.

Finally, the results of the analysis of the lattice data set with
$M_{\pi}=150$ MeV is reported in Fig. 8. The study of the vector form factor
suggest a branch point at time-like momenta but at $|p^{2}|$ that are smaller
than the branch point found for the $M_{\pi}=290$ MeV data. However, once
more, its precise location is difficult to determine within the statistical
precision of the simulation. On the other hand, both the differential
evolution and the simulated annealing methods point towards the presence of a
pole at $p^{2}=-0.114$ GeV2 with $Z=1.170$ for $N=2$ or $p^{2}=-0.094$ GeV2
with $Z=1.268$ for $N=4$ (differential evolution), or $p^{2}=-0.114$ GeV2 with
$Z=1.170$ for $N=2$ and $p^{2}=-0.117$ GeV2 with $Z=1.179$ for $N=3$
(simulated annealing), respectively. In what concerns the analysis of the
scalar form factor, the results suggest the presence of branch point at
similar values of time-like momenta as the found for the vector form factor
data. However, once more, its precise location is hard to define. This pole
occurs for $p^{2}=-0.153$ GeV2 with a $Z=0.468$ for $N=8$ and $p^{2}=-0.122$
GeV2 with $Z=0.368$ for $N=10$, according to the differential evolution
method, and at $p^{2}=-0.161$ GeV2 with $Z=0.463$ for $N=4$ and $p^{2}=-0.170$
GeV2 for $Z=0.522$ for $N=5$ when using the simulated annealing method.
Gathering these results together, one can claim a pole at time-like momenta
$p^{2}=-0.11$ GeV2 with a residuum $Z=1.17$ from the vector form factor data,
and a pole at $p^{2}=-0.15$ GeV2 for $Z=0.45$ from the scalar form factor
data. Once more, the residuum of the quark propagator at the pole is positive
defined.

Figure 9: Pion mass squared as a function of the pole quark mass as measured
from vector and scalar form factors. The dotted line refers to the physical
pion mass taken as $M_{\pi}=140$ MeV. To highlight the linear dependence of
the pion mass squared on the quark mass we include the solid lines that are
built using the average values of the quark masses for each of the pion mass
values and performing a linear regression (vector form factor in orange) or
going through the data points (scalar form factor in blue).

The results just described suggest also that there is a correlation between
the pion mass and the quark mass taken from the dominant pole, i.e. the pole
with the largest residuum, for real on-axis momenta. Defining the quark mass
as $m^{2}_{q}=-p^{2}$, where $p^{2}$ is the value of pole position for real
on-axis momenta, this correlation can be used to check one of the most
fundamental relations of QCD, namely the prediction based on the partial
conservation of the axial current, see e.g. Yndurain:2006amm , that gives

$M^{2}_{\pi}~{}\propto~{}m_{q}\ .$ (6)

The function $M^{2}_{\pi}(m_{q})$ is reported in Fig. 9 for the various
estimates of $m_{q}$ mentioned previously and it also includes the physical
value of the pion mass taken to be 140 MeV. For the various $m_{q}$ associated
with a given $M_{\pi}$ we take the pole mass for the smaller $N$ in the Padé
sequences. There are two key points that can be read from Fig. 9. The first
being that, despite the large errors, the data reported is compatible with the
predictions of PCAC given in Eq. (6). The other point being that the pole mass
for the scalar and vector form factors seems to be slightly different, with
the pole position of the scalar form part of the propagator appearing at
slightly higher values of $m_{q}$.

##  V Summary and Conclusions

The investigation of the poles and branch cuts using Padé approximants for the
Landau gauge pure Yang-Mills gluon and ghost lattice propagators and for the
full QCD quark propagator show significant differences between them. Indeed,
only for the gluon case the method clearly identifies poles at complex
momenta. The presence of complex poles associated with the gluon propagator
were already seen in Dudal:2010tf , where the compatibility of the predictions
of Gribov-Zwanziger type of actions, see Vandersickel:2012tz  and references
therein, with the propagators calculated with lattice QCD simulations was
investigated. In this sense, the observation of complex poles for the gluon
propagator favors the Gribov-Zwanziger type of approach to the quantization of
non-Abelian theories but it also rises questions about the definition of a
proper quantum theory for non-Abelian theories; see e.g. the discussion in
Hayashi:2021nnj ; Hayashi:2021jju  and references therein. The absence of
complex conjugate poles for the pure Yang-Millgs ghost and full QCD quark
propagators requires understanding the dynamics of the theory. If only one of
the fundamental QCD propagators has poles at complex momenta, it implies that
a delicate tunning has to take place that prevents complex poles in the ghost
and quark propagators. Looking at the Dyson-Schwinger equations for the quark
propagator and for the ghost propagtor, the mechanism responsible for the
absence of the complex poles for these two propagator should translate in some
type of constrains for the full ghost-gluon vertex and the full quark-gluon
vertex. Recall that the aforementioned Dyson-Schwinger equations involve only
these two vertices besides the propagators.

The known solutions of the Dyson-Schwinger equations for the quark propagator
suggest the presence of poles at complex momenta that the Padé method is
unable to identify. If this is a limitation of using Padé approximants to look
at the analytic structure of the propagators, or it is due to the statistical
precision of the simulations, or a limitation of the available analysis of the
Dyson-Schwinger gap equation remains to be investigated. Despite this
difference, the quark and the ghost propagators computed from these equations
are in good agreement with the outcome of the lattice simulations, suggesting
that the description of the fundamental QCD propagators with continuum
functional methods captures the essential of the dynamics of the theory.

The Padé analysis of the Landau gauge quark propagator identifies also a pole
for Minkowski momentum, i.e. real on-axis negative Euclidean momentum, with a
positive residuum. The naive interpretation of such a pole translates into a
free asymptotic single quark state. However, it is well known that positive
violation occurs for the quark propagator, i.e. the spectral density of the
quark is not positive defined, and the quark propagator can only be described
by a single pole in addition with other structures that, overall, prevent the
quarks to appear as free particles. Unfortunately, the method used herein does
not provide information on the structure beyond the single pole of the
fermionic propagator.

The quark propagator pole position in momentum space can be translated into an
effective quark mass. This effective mass falls is in the range of values
that, typically, are associated with the effective quark mass used in quark
models. As can be seen in Fig. 9, this effective quark mass is above 300 MeV.
Moreover, the same Fig. suggests that the pole position can differ for the
scalar and vector form factors, with the scalar form factors favouring
slightly large values for the effective mass. This result is somehow
unexpected and maybe an artefact of the method. Indeed, assuming that there is
a spectral representation for the quark propagator, then it can be written as
(ignoring the color part)

$S(p)=\int_{-\infty}^{+\infty}~{}dw~{}\frac{\rho(w)}{\not{p}-w+i\epsilon}\ ,$
(7)

where $\rho(w)$ is the quark spectral function Delbourgo:1977jc . From the
above expression it is difficult to accommodate a different pole structure for
the vector and scalar parts of the propagator. Alternatively, this could be an
indication of the absence of such a type of integral representation for the
propagator. One should also take into consideration the large uncertainties on
the estimation of the effective mass and that it is not clear on how to
estimate the errors for the effective mass.

The results of the Padé analysis also show a correlation between the effective
quark mass and the corresponding pion mass. Furthermore, it is shown in Fig. 9
that this correlation, measured by the curve $M^{2}_{\pi}(m_{q})$, is
compatible with the predictions of partial conservation of the axial current
for QCD.

The work described in the current manuscript and in Falcao:2020vyr  for the
pure gauge theory propagators shows that one can rely on Padé approximants
combined with global optimisation techniques to access the analytic structure
of the two point correlation functions. The accuracy on the results depends on
the size of the statistical ensemble of configurations and, certainly, having
access to the simulations with larger sets of gauge configurations will clear
the outcome of the analysis.

###### Acknowledgements.

The authors where supported by national funds from FCT – Fundação para a
Ciência e aTecnologia, I. P., within the projects UIDB/04564/2020,
UIDP/04564/2020. A.F.F. acknowledges the financial support via the Starting
Grant from Trond Mohn Foundation (BFS2018REK01) and the University of Bergen.

## References

  * (1) R. Alkofer and L. von Smekal, Phys. Rept. 353, 281 (2001) doi:10.1016/S0370-1573(01)00010-2 [arXiv:hep-ph/0007355 [hep-ph]].
  * (2) C. S. Fischer, J. Phys. G 32, R253-R291 (2006) doi:10.1088/0954-3899/32/8/R02 [arXiv:hep-ph/0605173 [hep-ph]].
  * (3) D. Binosi and J. Papavassiliou, Phys. Rept. 479, 1-152 (2009) doi:10.1016/j.physrep.2009.05.001 [arXiv:0909.2536 [hep-ph]].
  * (4) M. L. Perl, E. R. Lee and D. Loomba, Ann. Rev. Nucl. Part. Sci. 59, 47-65 (2009) doi:10.1146/annurev-nucl-121908-122035
  * (5) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, no.3, 030001 (2018) doi:10.1103/PhysRevD.98.030001
  * (6) J. Papavassiliou, [arXiv:2207.04977 [hep-ph]].
  * (7) D. B. Leinweber et al. [UKQCD], Phys. Rev. D 60, 094507 (1999) [erratum: Phys. Rev. D 61, 079901 (2000)] doi:10.1103/PhysRevD.60.094507 [arXiv:hep-lat/9811027 [hep-lat]].
  * (8) D. Becirevic, P. Boucaud, J. P. Leroy, J. Micheli, O. Pene, J. Rodriguez-Quintero and C. Roiesnel, Phys. Rev. D 60, 094509 (1999) doi:10.1103/PhysRevD.60.094509 [arXiv:hep-ph/9903364 [hep-ph]].
  * (9) D. Becirevic, P. Boucaud, J. P. Leroy, J. Micheli, O. Pene, J. Rodriguez-Quintero and C. Roiesnel, Phys. Rev. D 61, 114508 (2000) doi:10.1103/PhysRevD.61.114508 [arXiv:hep-ph/9910204 [hep-ph]].
  * (10) A. Cucchieri and T. Mendes, PoS LATTICE2007, 297 (2007) doi:10.22323/1.042.0297 [arXiv:0710.0412 [hep-lat]].
  * (11) I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker and A. Sternbeck, Phys. Lett. B 676, 69-73 (2009) doi:10.1016/j.physletb.2009.04.076 [arXiv:0901.0736 [hep-lat]].
  * (12) D. Dudal, O. Oliveira and N. Vandersickel, Phys. Rev. D 81, 074505 (2010) doi:10.1103/PhysRevD.81.074505 [arXiv:1002.2374 [hep-lat]].
  * (13) A. Cucchieri, D. Dudal, T. Mendes and N. Vandersickel, Phys. Rev. D 85, 094513 (2012) doi:10.1103/PhysRevD.85.094513 [arXiv:1111.2327 [hep-lat]].
  * (14) A. Maas, J. M. Pawlowski, L. von Smekal and D. Spielmann, Phys. Rev. D 85, 034037 (2012) doi:10.1103/PhysRevD.85.034037 [arXiv:1110.6340 [hep-lat]].
  * (15) D. Dudal, O. Oliveira and P. J. Silva, Annals Phys. 397, 351-364 (2018) doi:10.1016/j.aop.2018.08.019 [arXiv:1803.02281 [hep-lat]].
  * (16) S. W. Li, P. Lowdon, O. Oliveira and P. J. Silva, Phys. Lett. B 803, 135329 (2020) doi:10.1016/j.physletb.2020.135329 [arXiv:1907.10073 [hep-th]].
  * (17) G. T. R. Catumba, O. Oliveira and P. J. Silva, Phys. Rev. D 103, no.7, 074501 (2021) doi:10.1103/PhysRevD.103.074501 [arXiv:2101.04978 [hep-lat]].
  * (18) P. Boucaud, J. P. Leroy, A. Le Yaouanc, A. Y. Lokhov, J. Micheli, O. Pene, J. Rodriguez-Quintero and C. Roiesnel, JHEP 03, 076 (2007) doi:10.1088/1126-6708/2007/03/076 [arXiv:hep-ph/0702092 [hep-ph]].
  * (19) M. Q. Huber, R. Alkofer, C. S. Fischer and K. Schwenzer, Phys. Lett. B 659, 434-440 (2008) doi:10.1016/j.physletb.2007.10.073 [arXiv:0705.3809 [hep-ph]].
  * (20) D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel and H. Verschelde, Phys. Rev. D 78, 125012 (2008) doi:10.1103/PhysRevD.78.125012 [arXiv:0808.0893 [hep-th]].
  * (21) A. C. Aguilar, D. Binosi and J. Papavassiliou, Phys. Rev. D 78, 025010 (2008) doi:10.1103/PhysRevD.78.025010 [arXiv:0802.1870 [hep-ph]].
  * (22) J. Rodriguez-Quintero, JHEP 01, 105 (2011) doi:10.1007/JHEP01(2011)105 [arXiv:1005.4598 [hep-ph]].
  * (23) A. C. Aguilar, D. Binosi and J. Papavassiliou, Phys. Rev. D 81, 125025 (2010) doi:10.1103/PhysRevD.81.125025 [arXiv:1004.2011 [hep-ph]].
  * (24) S. Strauss, C. S. Fischer and C. Kellermann, Prog. Part. Nucl. Phys. 67, 239-244 (2012) doi:10.1016/j.ppnp.2011.12.025
  * (25) M. Peláez, M. Tissier and N. Wschebor, Phys. Rev. D 90, 065031 (2014) doi:10.1103/PhysRevD.90.065031 [arXiv:1407.2005 [hep-th]].
  * (26) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, J. Rodríguez-Quintero and S. Zafeiropoulos, Eur. Phys. J. C 80, no.2, 154 (2020) doi:10.1140/epjc/s10052-020-7741-0 [arXiv:1912.12086 [hep-ph]].
  * (27) U. Reinosa, J. Serreau, R. C. Terin and M. Tissier, SciPost Phys. 10, no.2, 035 (2021) doi:10.21468/SciPostPhys.10.2.035 [arXiv:2004.12413 [hep-th]].
  * (28) C. S. Fischer and M. Q. Huber, Phys. Rev. D 102, no.9, 094005 (2020) doi:10.1103/PhysRevD.102.094005 [arXiv:2007.11505 [hep-ph]].
  * (29) G. Eichmann, J. M. Pawlowski and J. M. Silva, Phys. Rev. D 104, no.11, 114016 (2021) doi:10.1103/PhysRevD.104.114016 [arXiv:2107.05352 [hep-ph]].
  * (30) M. Peláez, U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Rept. Prog. Phys. 84, no.12, 124202 (2021) doi:10.1088/1361-6633/ac36b8 [arXiv:2106.04526 [hep-th]].
  * (31) D. Dudal, D. M. van Egmond, U. Reinosa and D. Vercauteren, [arXiv:2206.06002 [hep-th]].
  * (32) L. von Smekal, R. Alkofer and A. Hauck, Phys. Rev. Lett. 79, 3591-3594 (1997) doi:10.1103/PhysRevLett.79.3591 [arXiv:hep-ph/9705242 [hep-ph]].
  * (33) P. Boucaud, J. P. Leroy, A. Le Yaouanc, A. Y. Lokhov, J. Micheli, O. Pene, J. Rodriguez-Quintero and C. Roiesnel, Eur. Phys. J. A 31, 750-753 (2007) doi:10.1140/epja/i2006-10295-1 [arXiv:hep-ph/0701114 [hep-ph]].
  * (34) R. Alkofer, W. Detmold, C. S. Fischer and P. Maris, Phys. Rev. D 70, 014014 (2004) doi:10.1103/PhysRevD.70.014014 [arXiv:hep-ph/0309077 [hep-ph]].
  * (35) R. Alkofer, C. S. Fischer, H. Reinhardt and L. von Smekal, Phys. Rev. D 68, 045003 (2003) doi:10.1103/PhysRevD.68.045003 [arXiv:hep-th/0304134 [hep-th]].
  * (36) A. G. Duarte, O. Oliveira and P. J. Silva, Phys. Rev. D 94, no.1, 014502 (2016) doi:10.1103/PhysRevD.94.014502 [arXiv:1605.00594 [hep-lat]].
  * (37) P. Boucaud, F. De Soto, J. Rodríguez-Quintero and S. Zafeiropoulos, Phys. Rev. D 96, no.9, 098501 (2017) doi:10.1103/PhysRevD.96.098501 [arXiv:1704.02053 [hep-lat]].
  * (38) A. G. Duarte, O. Oliveira and P. J. Silva, Phys. Rev. D 96, no.9, 098502 (2017) doi:10.1103/PhysRevD.96.098502 [arXiv:1704.02864 [hep-lat]].
  * (39) G. Krein, C. D. Roberts and A. G. Williams, Int. J. Mod. Phys. A 7, 5607-5624 (1992) doi:10.1142/S0217751X92002544
  * (40) P. O. Bowman, U. M. Heller and A. G. Williams, Phys. Rev. D 66, 014505 (2002) doi:10.1103/PhysRevD.66.014505 [arXiv:hep-lat/0203001 [hep-lat]].
  * (41) P. Boucaud, F. de Soto, J. P. Leroy, A. Le Yaouanc, J. Micheli, H. Moutarde, O. Pene and J. Rodriguez-Quintero, Phys. Lett. B 575, 256-267 (2003) doi:10.1016/j.physletb.2003.08.065 [arXiv:hep-lat/0307026 [hep-lat]].
  * (42) C. S. Fischer and M. R. Pennington, Phys. Rev. D 73, 034029 (2006) doi:10.1103/PhysRevD.73.034029 [arXiv:hep-ph/0512233 [hep-ph]].
  * (43) S. Furui and H. Nakajima, Phys. Rev. D 73, 094506 (2006) doi:10.1103/PhysRevD.73.094506 [arXiv:hep-lat/0602027 [hep-lat]].
  * (44) S. Furui and H. Nakajima, Phys. Rev. D 73, 074503 (2006) doi:10.1103/PhysRevD.73.074503
  * (45) A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti and J. Rodriguez-Quintero, Phys. Rev. D 86, 074512 (2012) doi:10.1103/PhysRevD.86.074512 [arXiv:1208.0795 [hep-ph]].
  * (46) E. L. Solis, C. S. R. Costa, V. V. Luiz and G. Krein, Few Body Syst. 60, no.3, 49 (2019) doi:10.1007/s00601-019-1517-9 [arXiv:1905.08710 [hep-ph]].
  * (47) J. R. Lessa, F. E. Serna, B. El-Bennich, A. Bashir and O. Oliveira, [arXiv:2202.12313 [hep-ph]].
  * (48) O. Oliveira, P. J. Silva, J. I. Skullerud and A. Sternbeck, Phys. Rev. D 99, no.9, 094506 (2019) doi:10.1103/PhysRevD.99.094506 [arXiv:1809.02541 [hep-lat]].
  * (49) P. Maris, Phys. Rev. D 50, 4189-4193 (1994) doi:10.1103/PhysRevD.50.4189
  * (50) P. Maris, Phys. Rev. D 52, 6087-6097 (1995) doi:10.1103/PhysRevD.52.6087 [arXiv:hep-ph/9508323 [hep-ph]].
  * (51) C. J. Burden, Phys. Rev. D 57, 276-286 (1998) doi:10.1103/PhysRevD.57.276 [arXiv:hep-ph/9702411 [hep-ph]].
  * (52) M. S. Bhagwat, M. A. Pichowsky, C. D. Roberts and P. C. Tandy, Phys. Rev. C 68, 015203 (2003) doi:10.1103/PhysRevC.68.015203 [arXiv:nucl-th/0304003 [nucl-th]].
  * (53) T. Frederico, D. C. Duarte, W. de Paula, E. Ydrefors, S. Jia and P. Maris, [arXiv:1905.00703 [hep-ph]].
  * (54) Y. Hayashi and K. I. Kondo, Phys. Rev. D 103, no.11, L111504 (2021) doi:10.1103/PhysRevD.103.L111504 [arXiv:2103.14322 [hep-th]].
  * (55) Y. Hayashi and K. I. Kondo, Phys. Rev. D 104, no.7, 074024 (2021) doi:10.1103/PhysRevD.104.074024 [arXiv:2105.07487 [hep-th]].
  * (56) A. F. Falcão, O. Oliveira and P. J. Silva, Phys. Rev. D 102, no.11, 114518 (2020) doi:10.1103/PhysRevD.102.114518 [arXiv:2008.02614 [hep-lat]].
  * (57) C. Pommerenke, J. Math. Anal. Appl. 41, 775 (1973)
  * (58) G. A. Baker, “Essentials of Padé approximants”, Academic Press, U.S.A., (1975)
  * (59) G. A. Baker and P. Graves-Morris, “Encyclopedia of mathematics and its applications”, Cambridge University Press, Cambridge, U.K., (1996)
  * (60) C. Bender and S. Orszag, “Advanced mathematical methods for scientists and engineers I: asymptotic methods and perturbation theory”, Springer, New York, U.S.A., (1999)
  * (61) J. Sanz-Cillero, [arXiv:1002.3512 [hep-ph]].
  * (62) P. Masjuan Queralt, [arXiv:1005.5683 [hep-ph]].
  * (63) https://www.wolfram.com/mathematica/
  * (64) S. W. Li, P. Lowdon, O. Oliveira and P. J. Silva, Phys. Lett. B 823, 136753 (2021) doi:10.1016/j.physletb.2021.136753 [arXiv:2109.10942 [hep-th]].
  * (65) D. Binosi and R. A. Tripolt, Phys. Lett. B 801, 135171 (2020) doi:10.1016/j.physletb.2019.135171 [arXiv:1904.08172 [hep-ph]].
  * (66) D. Zwanziger, Nucl. Phys. B 323, 513-544 (1989) doi:10.1016/0550-3213(89)90122-3
  * (67) D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel and H. Verschelde, Phys. Rev. D 78, 065047 (2008) doi:10.1103/PhysRevD.78.065047 [arXiv:0806.4348 [hep-th]].
  * (68) D. Dudal, S. P. Sorella, N. Vandersickel and H. Verschelde, Phys. Rev. D 77, 071501 (2008) doi:10.1103/PhysRevD.77.071501 [arXiv:0711.4496 [hep-th]].
  * (69) O. Oliveira and P. J. Silva, Phys. Rev. D 86, 114513 (2012) doi:10.1103/PhysRevD.86.114513 [arXiv:1207.3029 [hep-lat]].
  * (70) A. Cucchieri, D. Dudal, T. Mendes and N. Vandersickel, Phys. Rev. D 93, no.9, 094513 (2016) doi:10.1103/PhysRevD.93.094513 [arXiv:1602.01646 [hep-lat]].
  * (71) F. J. Yndurain, “The Theory of Quark and Gluon Interactions,” doi:10.1007/3-540-33210-3
  * (72) N. Vandersickel and D. Zwanziger, Phys. Rept. 520, 175-251 (2012) doi:10.1016/j.physrep.2012.07.003 [arXiv:1202.1491 [hep-th]].
  * (73) A. C. Aguilar, D. Binosi, C. T. Figueiredo and J. Papavassiliou, Eur. Phys. J. C 78, no.3, 181 (2018) doi:10.1140/epjc/s10052-018-5679-2 [arXiv:1712.06926 [hep-ph]].
  * (74) R. Delbourgo and P. C. West, J. Phys. A 10, 1049 (1977) doi:10.1088/0305-4470/10/6/024