Spectral functions of confined particles

We determine the gluon and ghost spectral functions along with the analytic structure of the associated propagators from numerical data describing gauge correlators at space-like momenta obtained by either solving the Dyson-Schwinger equations or through lattice simulations. Our novel reconstruction technique shows the expected branch cut for the gluon and the ghost propagator, which, in the gluon case, is supplemented with a pair of complex conjugate poles. Possible implications of the existence of these poles are briefly addressed.


Introduction
As a consequence of the self-interactions of gauge fields, SU(N C ) non-Abelian gauge theories are asymptotically free, i.e., the coupling constant which controls the strength of their interactions decreases as the momentum scale increases [1,2]. This characteristic behavior persists even when matter particles are added, provided that their number N f is less than a critical value (11N C /2 at leading order, which is further reduced by non-perturbative effects, see, e.g. [3]). This is indeed the case for Quantum Chromo-Dynamics (QCD), the theory thought to describe the strongly interacting sector of the Standard Model, in which one has N C = 3 (resulting in 8 adjoint gluons) and N f = 6 (corresponding to 6 fundamental quark flavors). The converse is also true: the QCD interaction strength grows with the separation between gluons and quarks. The theory is confined: colored states, like isolated gluons and/or quarks, do not appear in its spectrum.
Knowledge on the full propagator is often equalled with having solved the underlying theory. It is the central object in the calculation of equilibrium as well as dynamical observables with applications ranging from the determination of the hadron and glueball spectrum to the calculation of transport coefficients of the quark gluon plasma. The determination of the full propagator and its analytic structure in the complex plane is therefore indispensable.
But how does the analytic structure of these correlators look like? Consider for simplicity the pure glue theory in which the only dynamical scale is Λ QCD . Then, perturbation theory implies that, whenever the condition q 2 ≡ q 2 /Λ 2 QCD 1 is met, one has at leading order [33] ∆(q) ∼ Z gl /q 2 (log q 2 ) γ and D(q) ∼ Z gh /q 2 (log q 2 ) δ where Z gl and Z gh are dimensionless renormalization constants and γ and δ are the (one-loop) anomalous dimensions (with γ = 13/22 and δ = 9/44 for N f = 0, which will be the case from now on). Accounting for the perturbative logarithmic running in the ultraviolet (UV) momentum region, requires the presence of a branch cut structure for real and time-like q 2 momenta. The presence of this cut, without any additional singularities in the complex plane, would in turn allow for the usual Källén-Lehmann representation by which the full propagator is expressed as an integral over the (free) propagator with a spectral density depending on the frequency ω.
In fact, owing to the axioms of local quantum field theory [34], any 2-point correlation function must be an analytic function in the cut complex q 2 -plane with singularities along the time-like real axis only [35]; for any other singularity structure to be present, e.g., (simple) poles, one or more of these axioms, typically local space-like commutativity, i.e., causality, needs to be relaxedAs has been discussed in [36][37][38], however, although causality might be violated at the level of the 2-point functions in such cases, the corresponding S -matrix remains both causal and unitary. One might also argue [15] that strict local- . Middle: the propagator curve plotted obtained from the smallest norm spectral function. Right: histogram of all poles of the SPM propagators that lie in the left half-plane, thus omitting Froissart doublets in the right half-plane; the counting frequency is defined as N/N max , where N is the number of poles in a given bin and N max the maximum number of poles per bin found in the histogram; the branch cut at q 2 ≤ 0 and the poles at q 2 1 = −1 ± i are correctly identified. The poles used in (1) are only those with residues above a certain threshold, see text for details. ity may be too strong a requirement for physical theories displaying confinement, in view of non-local complications like, e.g., the Gribov-Singer ambiguity [39,40]. Furthermore, several phenomenological models predict the appearance of complex simple poles (which are bound to materialize in conjugate pairs, for the propagator is real at space-like momenta) in the gluon propagator [12,41,42], see also [39,43,44]. Further evidence for their existence has been found in the quark sector where the Nakanishi representation [45][46][47] perfectly describes solutions of the quark gap equation [48] using the most sophisticated quark-gluon kernel available [49,50] (see also [51,52]).
Thus, to ascertain whether or not such complex conjugate poles appear in the analytic structure of QCD propagators, a method that assumes no (or at least as minimal as possible) prior knowledge on such structures is needed. There are several numerical continuation methods available in the literature that aim at obtaining the best possible reconstruction of spectral functions [53][54][55][56][57][58][59][60][61][62][63][64]; however, as all these reconstruction techniques rely on the standard Källen-Lehmann-type spectral representation, they ignore the possible presence of additional structures in the complex plane. The purpose of this letter is to describe a novel numerical analytic continuation tool allowing for this possibility, and apply it to the available continuum and lattice data for the gluon and ghost propagators in order to obtain their respective spectral functions together with the corresponding analytic structure.

Generalized spectral representation and reconstruction method
In the presence of n complex conjugate pairs of simple poles located at q i , q * i (i = 1, .., n) and with residue R i , R * i , the usual Källén-Lehman spectral representation gets generalized to where q refers now to the energy-component and ρ D is the (nonpositive definite) spectral density defined at real frequencies ω, with → 0 + [65,66]. As D(q) ∼ 1/q 2 in the IR, for the ghost particle we will study the spectral function associated to its dressing function ρ F rather than ρ D , the two being related by [64]. We also note that assuming q 2 ∆(q) → 0 for q 2 → ∞, as required by perturbation theory, implies a modified version of the usual Oehme-Zimmermann superconvergence (OZS) relation [67,68], which now becomes The proposed reconstruction algorithm is a variant of the Schlessinger Point Method (SPM) which is based on a rationalfraction representation similar to Padé approximants [69]. Given a set of N input points (q i , y i = D(q i )) one first constructs a continued fraction representation of the propagator: where the coefficients a i are recursively evaluated such that the function D(x) acts as an interpolation through all the points, i.e., D(q i ) = y i , i = 1, 2, . . . , N (see [69][70][71] for details). Once the function D(q) is determined, its analytic continuation is defined as the meromorphic function obtained by allowing the originally real q to take on complex values; the corresponding spectral function, ρ D , can be then evaluated using Eq. (2). The quality of the reconstructed spectral function can be assessed by inserting ρ D in Eq. (1) and comparing the resulting propagator D ρ with the original one by means of a suitable norm, the simplest one being Notice that the OZS relation is not used to further constrain the spectral function since we focus on intermediate energies and the reconstructed spectral function by construction behaves asymptotically as 1/ω, giving rise to a divergent integral. In the following we will be dealing with data sets comprising M > N values for the propagators D in the space-like Euclidean momentum region, from which we want to reconstruct the associated spectral functions. We are therefore confronted with the problem of which N points should be chosen as our input. In fact, if the data were exact, i.e., without any statistical error, any subset chosen would give, in general, the same result. In realistic situations, however, data do have errors, and some subsets of input points will give results that are closer to the correct one (i.e., the one that would be obtained for exact data) than others. In order to identify such 'optimal' subsets we have devised an algorithm comprising the following steps: (i) Select a partition of N random points chosen from the complete set of M points (q i , y i ); (ii) Use Eq. (3) to get the corresponding D(q i ); (iii) Use Eq. (2) to calculate the spectral function; if there are known constraints, e.g., on its asymptotic behavior and/or its positiveness (or lack thereof), this information can be used here to either accept or reject the spectral function; (iv) Reconstruct the propagator D ρ from Eq. (1), taking only complex conjugate poles with a residue above a certain threshold into account, and evaluate its norm (4); (v) Repeat steps (i) through (iv) for L times and identify the input point (q j , y j ) that was used most often among the subset of K functions with the smallest norm; (vi) Repeat steps (i) through (v) N − 1 times, keeping fixed the identified points (q j , y j ) until all N optimal points have been selected.
As described in (iv), only poles with residues above a certain threshold are taken into account when evaluating Eq. (1). This threshold is chosen such as to discard non-physical Froissartdoublet poles while at the same time keeping candidates for physically meaningful poles. In the present work, we use a threshold in the range between 0.1 and 2. Let us now briefly discuss the values chosen for the different parameters. To begin with we note that, if D in Eq. (1) is expressed as a rational fraction, D = P/Q, P and Q will be polynomials of order (N − 1)/2 (P and Q) for an odd number of input points, and of order N/2 − 1 (P) and N/2 (Q) for an even number of input points. Even though knowing that the propagator vanishes in the UV leads to the natural choice of N even, one cannot hope in any case to numerically reproduce its asymptotic behavior correctly based on noisy input data; and, in fact, in the intermediate energy range considered here, we checked that our results are stable irrespectively of N being even or odd. What matters is that there seems to be always an optimal number of input points, call it N * : for fewer input points we found that the reconstruction depends strongly on the number of points, whereas a larger number gives rise to numerical instabilities due to a loss of accuracy when calculating the coefficients of the continued fraction (2) (the C++ code developed uses the GMP library for arbitrary precision arithmetic on floating point numbers). For all datasets we have found that choosing around 50 input points always works, and therfore set N * = 50. For the remaining parameters we have set L = 4000 and K = 200 which offered a good compromise between execution speed and results precision. M varies among the different data sets from a minimum of 69 (DS gluon case) to a maximum of 148 (lattice ghost case); finally we effectively set appearing in Eq. (2) equal to zero.
As step (i) of the algorithm is intrinsically random, every run of the algorithm will result, in general, in a different set of optimal points. Typically, we had 500-1000 full runs per case (performed on the Goethe-HLR cluster), out of which we selected only the ones with a norm below a suitable threshold, set at 1.5 times the minimum value reached. All solutions are then resampled in order to construct the mean and the standard deviation, thus furnishing respectively the best curve and the one and two σ confidence levels (gray error bands in the figures).
We proceed with a model application of the algorithm in order to demonstrate its capabilities. For this purpose, we generate mock propagator data consisting of 100 points (q i , y i ) in the range q i ∈ [0.01, 10] GeV, using Eq. (1) with a Breit-Wigner spectral function 2Γω/[π((ω 2 − Γ 2 − M 2 ) 2 + 4Γ 2 ω 2 )] with M = 4Γ = 1 GeV and supplemented with a pair of complex-conjugate poles located at q 2 1 = −1 ± i GeV 2 with residue R 1 = 1. Finally, we add a statistical error to the input data through y i → y i (1 + εr i ) with ε = 10 −3 and r i a random number drawn from a normal distribution with zero mean and unit standard deviation. We note that only the central values of the data points are used for the reconstruction. The uncertainties are therefore taken into account effectively, as encoded in the scattering of the central values. The error range can in principle also be taken into account explicitly by generating data points within the given error bars. As can be seen in Fig. 1, the method determines the correct ρ D and the corresponding propagator D ρ , also accurately identifying the location of the poles in the complex q 2 -plane, see also Table 1. The presence of a branch cut, to be expected on analytic grounds due to a square root term appearing when integrating the Breit-Wigner spectral function, is also visible as a sequence of poles on the negative real axis. Indeed, the reconstruction of the meromorphic structure of the function is possibly the most stable result that the algorithm gives: increasing the error rate ε results in a loss of precision in locating peak position and height of ρ D , but one will still get a reliable estimate on the cut and poles positions.

Results
Our results for the gluon spectral function 1 ρ ∆ are shown in Fig. 2, in which the upper panels correspond to input points 1 To ease the comparison between different results, we normalize both the DS and lattice data so that ∆(0) = 1 [GeV −2 ] and F(0)=1. This should not be confused with renormalization, e.g., it is known that it is not possible to renormalize F to be unity at q = 0 without violating symmetry identities [72]. obtained from the numerical solution of the Dyson-Schwinger (DS) system for the gluon and ghost two-point functions [16], whereas the lower panels use as input the 64 4 at β = 6.0 lattice results for SU(3) Yang-Mills theory of [8]. Qualitatively the spectral functions behave in the same way, approaching large and small frequencies from the negative region as required [73], and displaying the same peak sequence; in addition, the negative contributions to the spectral function confirm the confined nature of the gluon, removing it from the theory asymptotic spectrum. The reconstructed propagators obtained from Eq. (1) reproduce the input data well (with residuals at the 10 −3 level, see insets in the middle panels). Perhaps, the most striking resemblance between the two cases can be found in the reconstructed analytic q 2 structure, featuring a branch cut singularity at the ghost threshold Re q 2 ≤ 0 and a single pair of complex conjugate poles located at q 2 1 ≈ −0.2 ± 0.3i GeV 2 (DS) and q 2 1 ≈ −0.3 ± 0.5i GeV 2 (lattice), see Table 1 again. The presence of these poles also makes clear the quantitative difference between the two spectral functions: since the residue is larger in the lattice case with respect to the DS one, more weight for the reconstruction of the propagator is carried by the poles in the former case. Notice that these results are qualitatively different from the ones reported in [16] (DS) and [64] (lattice), where the only analytic structure found was a branch cut. Knowledge of ρ ∆ allows for defining a renormalization group invariant effective gluon mass m gl . Even though this notion has proven useful, see, e.g., [9,16,[74][75][76][77][78][79][80][81][82], it should be clear that m gl is not an observable (and the gluon not a regular massive particle); rather, m gl coincides with the scale at which positivity violation appears. From the zero crossing of the spectral function (left panels insets if Fig. 2) we can estimate m gl ≈ 0.4 ± 0.1 GeV (DS) and m gl ≈ 0.6±0.2 GeV (lattice). These values are compatible with available theoretical and phenomenological estimates of the RG effective gluon mass [9,16,[74][75][76][77][78][79][80][81], and, in particular, the estimate m gl ≈ 0.5 GeV obtained when constructing the process-independent QCD effective charge [72,83] analogue of the quantum electrodynamics Gell-Mann-Low effective coupling [82].
We note that a qualitatively similar behavior is observed when analyzing data on the gluon propagator obtained in [21]. This data corresponds to a Landau gauge gluon propagator of the so-called scaling type, i.e., a type of solution of the DS system that is characterized by the IR power laws ∆(q) ∼ (q 2 ) 2κ−1 and F(q) ∼ (q 2 ) −κ with κ ≈ 0.58 [84]. While in [73] a modified Breit-Wigner Ansatz was used to reconstruct the gluon spectral function, which by definition excludes the presence of poles in the q 2 plane, applying the SP method yields a clear signal for complex conjugate poles at q 2 1 ≈ −0.1 ± 1.3i GeV 2 as well as a spectral function qualitatively similar to Fig. 2.
Turning to the ghost sector, the spectral function for the dressing function ρ F is shown in Fig. 3, with the upper (lower) panels corresponding to the DS (lattice) data. We find a qualitatively similar behavior to the gluon case, albeit with larger uncertainties. Notice also that the IR asymptotics obtained from lattice data looks opposite to the one coming from the DS data; however, the latter should be considered more reliable, due to the smaller uncertainties. The analytic structure shows the presence of only a branch cut (barring some spurious poles in the lattice case due to their lower precision), starting at zero as was already the case for the gluon. This has to be expected, due to the presence of the ghost gluon vertex and the fact that the gluon, being no ordinary massive particle, carries spectral weight down to arbitrarily small energies, see Fig. 2. Similar results (shown in Fig. 4) are obtained when analyzing the gluon and ghost two-point functions obtained from lattice configurations involving two light (twisted mass) degenerate quark flavours and two heavy ones [7]. The inclusion of dynamical quarks has been known to lead to a lower saturation point of the gluon propagator (see, e.g., Fig. 1 in [7]), which can be in turn interpreted as the gluon becoming "more massive" (an effect that has to be expected on theoretical grounds [85,86]); on the other hand, due to the absence of a direct coupling between ghosts and quarks, the ghost is practically unaffected by the presence of the latter particles. Indeed, the gluon spectral function (see Fig. 4 again) reveals that the renormalization group invariant gluon mass increases to m gl ≈ 0.75 ± 0.05 GeV; in addition, as was the case for the quenched data, one finds that the gluon spectral function is characterized by the presence of a single pair of complex conjugated poles located at q 2 1 ≈ −0.6 ± 0.5i GeV 2 (see Table 1), whereas the ghost shows only the logarithmic cut.

Discussion and conclusions
The strength of the signal in our numerical procedure coupled to the fact that the same pole singularities appear when analyzing both the DS and lattice data, suggest that their presence is a characteristic of the gluon 2-point function rather than an artifact. This calls in turn for a thorough reassessment of DS truncation schemes as well as gauge-fixing procedures for lattice-regularized simulations (it is known in the former case that changing the vertex in the quark gap equation can significantly alter the singularity structure of the solutions, see e.g., [87]). In this respect, however, we have considered the gluon propagator data of [25], which were obtained from lattice configurations gauge fixed in a covariant gauge at non-zero values of the gauge fixing parameter ξ. The results for the case ξ = 0.5, shown in Fig. 5, are of exactly the same type of the ones found in the Landau gauge case: a pair of complex conju-gate poles and a logarithmic cut 2 .
Also, in Ref. [88] it has been shown that the quenched Landau gauge lattice data of [8] can be described (in the Re(q 2 ) > 0 region) in terms of simple and double poles located on the real negative q 2 -axis. These poles (which, incidentally, would not break locality) would in turn correspond to unconventional singular terms (derivatives of the δ-function) present in the spectral representations of unphysical correlators, and which are permitted due to the indefinite metric characterizing covariant gauge theories (see again [88] and references therein). However, the analysis of mock data generated according to these predictions, shows that our method does correctly reconstruct such poles on the real negative q 2 -axis; and, in particular, that it does not mistakenly reconstruct them as complex conjugate pairs, even when noise is added. Therefore, the lack of any unambiguous signal for such structures in all the cases analyzed in this work (not only the quenched Landau gauge data of [8]) strongly suggest that they are not part of a veracious spectral representation describing confined particles.
Overall, these results seem to suggest that the presence of a complex conjugate poles pair is a characteristic feature of the gluon 2-point function in covariant gauges. If confirmed, this fact would then have far-reaching consequences. At the theoretical level one would be tempted to conclude that Yang-Mills theories break locality non-perturbatively; whereas, at the phenomenological level, it would provide a new microscopic foundation for building better models of the gluon and ghost propagators leading to the development of quark-gluon kernels with the potential of describing QCD bound-states, form factors, parton distribution amplitudes and functions, at the degree of precision required by upcoming experimental efforts. DS and lattice propagator data. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project number 315477589 TRR 211. The work of R.-A. T. is also supported by the BMBF under grant No. 05P18RFFCA.