Pad\'e approximation and glueball mass estimates in 3d and 4d with N_c = 2,3 colors

A Pad\'e approximation approach, rooted in an infrared moment technique, is employed to provide mass estimates for various glueball states in pure gauge theories. The main input in this analysis are theoretically well-motivated fits to lattice gluon propagator data, which are by now available for both SU(2) and SU(3) in 3 and 4 space-time dimensions. We construct appropriate gauge invariant and Lorentz covariant operators in the (pseudo)scalar and (pseudo)tensor sector. Our estimates compare reasonably well with a variety of lattice sources directly aimed at extracting glueball masses.


Introduction
Although confinement is a well accepted phenomenon in pure gauge theories [1], the extraction of the observable degrees of freedom, which ought to be glueballs, is a challenging task. Several theoretical methods 1 have been tested, to name a few: qualitative studies [2], effective Hamiltonian methods [3,4,5,6], AdS/CFT inspired tools [7,8], lattice simulations [9,10,11,12], functional approaches [13,14], Regge trajectory analyses [15], sum rules analyses [16], etc. We refer to [17] for a recent review on the subject. Also, from the experimental side, the status of glueballs is at the best inconclusive as they are hard to detect, partially due to their mixing with other states, see [18] for more details.
In this current note, we will apply a method developed in [19]. The main purpose is to benefit from high precision lattice computations of the gluon propagator in certain preferential gauges, in particular the Landau gauge [20]. As these correlation functions carry essential nonperturbative information on the gluon dynamics, it seems natural to benefit from these data. As we are interested in continuum computations, we need functional forms for e.g. the gluon propagator, as we plan to study the correlation functions of bound state and thus of composite operators. As it will become clear, we will probe the analyticity properties of the latter propagators, which is not an easy task if the input gluon propagator is a complicated function. A recent numerical approach to derive the spectral density of a (tree level) bound state propagator given an a priori analytical prescription for the input constituent propagator can be found in [21,22], thereby confirming the analytical results of [23]. A far more appealing approach would be to only use the gluon lattice data, which is however still in its infancy given all the difficulties to extend the lattice data from the Euclidean region p 2 ≥ 0 to the complex p 2 plane [24].
An effective action formalism to deal with the Gribov issue was worked out in a series of papers by Gribov & Zwanziger [29,30,31,32], see also the recent work [33]. Basically, the domain of integration of the gauge fields in the Euclidean functional integral is further constrained to the first Gribov region Ω, whose boundary is the Gribov horizon, where the Faddeev-Popov operator attains the first vanishing eigenvalue [34]. In recent years, we included into the original derivation of [29,30,31,32], the dynamical effects of dimension two condensates [35,36,37,38], resulting in what is nowadays called the Refined Gribov-Zwanziger (RGZ) action.
Let us also notice that the relevance of dimension two condensates for certain nonperturbative effects in gauge theories was already realized in e.g. [39,40,41,42,43]. The analytic form of the tree level gluon propagator obtained in the Refined Gribov-Zwanziger theory is, in general, given by where m 2 , resp. M 2 and ρ (with c.c. ρ † ) are corresponding to the d = 2 gluon condensate 〈A 2 〉, resp other d = 2 condensates constructed from the additional fields present in the Gribov-Zwanziger action, let us refer the interested reader to the paper [38] for details. The quantity λ 4 is directly related to the restriction to the Gribov region Ω. If it happens that ρ is real, then the RGZ propagator can be rewritten into the form The d = 2 mass parameters M 2 1 , M 2 2 and M 2 3 are recombinations of m 2 , M 2 , ρ and λ 4 , whereby M 2 and ρ appear in a fixed combination [38]. We feel it is important to point out here that the form of the propagator (1) is in general predicted by the RGZ formalism. The mass parameter λ 4 is present because of the restriction to the Gribov region Ω, the others are related to a stabilization of the vacuum by means of d = 2 condensates. Thus this tree level analytic form is dictated by the underlying RGZ dynamics. In principle, these vacuum expectation values are determined by self-consistent gap equations, that is by minimization of the effective potential (vacuum energy). Such program was carried out first in [36,38], albeit in a one-loop approximation, giving a qualitative but not always a superb quantitative agreement with the data. This is acceptable, since it is evidently not a straightforward task to compute the effective potential to arbitrary high order. The key features of the formalism are clear at lowest order already. Though, for the current purposes we also require a quantitative estimate for the condensates, as these are exactly the mass parameters that will fuel our eventual glueball mass estimates. We will therefore make use of the gluon lattice data to attribute values to the condensates, in particular have the functions (1) and (2) been used as fitting proposals in the papers [25,28] on which we shall thus rely.
It is worthwhile to remember that other functional approaches exist for the study of Yang-Mills Green functions in the Landau gauge, in particular the Schwinger-Dyson equations that corresponds to the quantum equations of motion. Also these equations are impossible to solve exactly, but can be replaced after certain approximations by (numerically) solvable equations. These numerical propagators can then also be fitted with functions as in [44,45]. One can even directly attempt to construct numerical estimates for the gluon spectral function based on either the Schwinger-Dyson equations [46] or "inversion" of the lattice data [47], but in none of the aforementioned cases a closed analytical expression can be derived, one is always reduced to fitting the numerical result with some a priori completely free to choose function. This is different from the RGZ (loop expansion) approach, where the functional form are closed analytical expressions (albeit with the condensates' values fixed via the lattice data).
In Sect. 2, we will first describe which composite operators we need to describe specific glueball states and we will construct the spectral density of the corresponding two-point correlation function 〈 (p ) (−p )〉 in a first order approximation, also known as the Born approximation [21]. We furthermore list the values, as reported and discussed in other works, for the RGZ parameters appearing in either (1) or (2). In Sect. 3, we spend a few words on the infrared moments technique and how it can be used to get a first rough estimate of glueball masses. We end with a discussion in Sect. 4.

2
The gluon propagator input, the glueball operators and associated spectral densities

The RGZ gluon propagator
Let us first give some numbers we will need later on. For the SU (3) case, we rely on the RGZ fitting parameters of [25] 2 , where the expression (2) seems to be singled out by the data 3 : With these values, the propagator (2) displays two c.c. poles, namely at in appropriate GeV units. To obtain the location of the poles, we always use the central value of the fitting estimates as in (3). The one standard deviation errors on those are taken from the original papers and shown between parentheses.
This observation of c.c. poles is what lies at the heart of the i -particles setup introduced in [48]: the RGZ gluon propagator can be expressed in terms of a pair of "complex" particles with c.c. masses. Clearly, such degrees of freedom are not physical, hence there is no direct observable information in the Landau gauge gluon propagator.
As discussed in [28] for SU (2) lattice data 4 , the corresponding fitting parameters are given by leading to the roots (GeV units) [28, Table IV] For SU (2), there are also data available in d = 3. Interestingly, there, the full RGZ propagator (2) is needed to adequately describe the lattice data. For convenience, we write it in the following form (GeV units) The d = 3 gluon propagator thus displays next to two c.c. roots also a relatively small real root. Though, since its residue is negative, it neither describes a physical degree of freedom.

(Pseudo)scalar and (pseudo)tensor glueball operator in d = 4
As we wish to work in the context of local Lorentz-invariant 5 quantum field theory, we are looking for gauge invariant/Lorentz covariant operators, constructed in such a way that they described states with specific J PC quantum numbers. This occasionally necessitates the introduction of projection operators onto the desired subspace. For the scalar operator ( J PC = 0 ++ ), we can take 0++ = F 2 µν , while for the pseudoscalar state (J PC = 0 −+ ) the corresponding operator is given by 0 In the (pseudo)tensor sector, a little more effort is required. As the equations of motion derived from the (RGZ) action are somewhat different from those of the usual Yang-Mills action [31,36], the standard energy momentum tensor, αβ does not qualify anymore as it is not necessarily conserved 6 . In [19,49], we proposed 2. Extrapolated to infinite volume with β = 6.2.
3. We will always omit the necessary global renormalization factor as it will play no role in the current paper. 4. Obtained with V = 128 4 as the largest volume and with β = 2.2. 5. In practice, we employ an Euclidean space-time. 6. This does not mean the (RGZ) theory has no conserved energy-momentum tensor, it simply differs from t µν .
for the tensor (J PC = 2 ++ ) glueball; where P µν = ∂ 2 δ µν −∂ µ ∂ ν is the transverse projector. The foregoing operator 2++ µν is, irrespective of the equations of motion, conserved, symmetric and traceless and has the upshot to be directly proportional to t µν if it happens that ∂ µ t µν = 0.
The symmetric and conserved operator can be easily made traceless, without compromising its symmetry properties as well as its conservation, namely by passing to an operator suited to describe the J PC = 2 −+ state.

Scalar and tensor glueball operator in d = 3
In d = 3, the dual of F µν is an axial vector. We have not been able to write down in a Lorentz covariant notation a simple local operator that would describe the analogue of pseudoscalar/-tensor as in d = 4. Therefore, for the present work, we will limit ourselves to the scalar/tensor case where the operators 0++ and 2++ µν can be immediately employed.

2.4
Derivation of the spectral densities at lowest order Using the i -particle representation of the RGZ propagator [19], the necessary spectral representations of the afore described glueball correlators can be derived using the tools of [23]. We recall here that a physical particle propagator, e.g. 〈 0++ (p ) 0++ (−p )〉 ≡ (p 2 ), should be consistent with a Källén-Lehmann integral form, which in Euclidean conventions reads with ρ(t ) t ≥τ0 ≥ 0. Making p 2 a complex variable, eq. (11) describes an everywhere analytic function, with the exception of a branch cut along the negative real axis. Using Cauchy's basic theorem, the density ρ(t ) is proportional to the discontinuity along the axis, which due to the optical theorem means that it must be positive, at least when we are talking about physical observables. E.g. for a confined gluon this does not need to be the case. Now, given that the i -particles come in complex pairs and that we shall consider the single bubble approximation 7 to (p 2 ), we need to carefully consider which diagrammatic subsets of the full correlator can correspond to a real degree of freedom, in the current approximation at least. As studied into great detail in [48], if the 2 internal propagator lines of the bubble corresponds to a pair of i -particles with c.c. masses, and only then, the resulting contribution is consistent with the representation (11). A case by case check is then all that is needed to verify the positivity of the density ρ(t ). For the moment, a full-fledged approach that can consistently remove the unphysical pieces (from combining particle propagators not containing c.c. mass pairs for example) order by order is not yet available. We will adopt the working assumption that for now, we can stick with only retaining the physical contributions to the respective correlators. Recalling that, upon considering particles with masses squared m 2 ≥ 0 and M 2 ≥ 0, their one-loop composite two-point correlation function will develop a branch point at p 2 = −(m + M ) 2 , see e.g. [50], it is clear that complex valued branch points are to be expected when propagator with complex masses are combined. This does not occur when m 2 and M 2 are c.c. , as discussed in [23,48].
So, using the complex mass Cutkosky rules (see [23] for a discussion), we can derive for each of the above operators the physical piece of the glueball correlators. Notice that in d = 3, also the 2 Yukawa propagators with the wrong sign, cf. (7), will contribute to the physical piece of the bound state correlator, since the 2 minus signs will combine into a physical, positive signed, piece. We will refrain from writing down the tedious calculations here, but immediately list the final spectral densities we will continue to work with. 7. All considered operators have quadratic pieces in the gluon field, hence the lowest order contribution to the considered bound state propagators will always be a bubble diagram.
In d = 4, we found, up to irrelevant global prefactors, The Heaviside step function (x ) implements the threshold which is in all cases given by the expression τ 0 = 2 µ 2 + µ 4 + 2θ 4 . It remains to check under which conditions the above spectral functions are positive. For (12) and (13), it is easily checked that there are no real roots for t > τ 0 and that the functions are positive for that interval given that always θ 2 ≥ 0, µ 2 ≥ 0. In the case of (15), the only real root for t > τ 0 is given by t = 3µ 2 + 3 3µ 4 − 16θ 4 , at least when 3µ 4 ≥ 16θ 4 . It is then easily checked that this root is located below τ 0 if µ 4 ≤ 98 15 θ 4 . If the latter condition is fulfilled, then (15) remains positive over the whole of the relevant domain. Using the numbers quoted in eq. (3) for SU (3) and eq. (5) for SU (2), the foregoing conditions are clearly fulfilled. For the tensor case (14), we had to check numerically that the spectral density is positive for t > τ 0 .
In d = 3, we found, up to irrelevant global prefactors, In both cases, no real roots are to be reported over the t -domain of interest for µ 2 ≥ 0, θ 2 ≥ 0, hence the d = 3 spectral densities are found positive.

Setup of the moment problem
So far, we have derived the spectral densities of a set of glueball correlators at lowest (one bubble) order. Clearly, none of the functions (12), (13), (14), (15) or (16), (17) displays a pole on the negative real axis, something that would correspond to a massive physical particle. This is no surprise given the approximation, usually a dynamical pole will only emerge if some kind of resummation is performed. Though, given the complexity of the gluon interactions, to our knowledge, a self-consistent approximation scheme that would allow to resum to some extent higher order bubble graphs, let stand alone construct the spectral properties of such graphs 8 , has not yet been achieved so far.
Therefore, we will adopt a different strategy here, based on a suitable moment problem. Moment problems are not new in particle physics, see [16] for similar approaches using the OPE/spectral methods/sum rules. The difference with what we will explain is that our moments will be IR based, in contrast with ruling approaches. Let us first provide a short survey of the IR moment problem as originally introduced in [23]. We reconsider (p 2 ) as given in eq. (11), and perform the substitution t = 1 s , so that 8. For the quark dynamics, a popular approximation scheme is that of Nambu-Jona-Lasinio-like (NJL) interactions based on integrating out the gluons, which allows resummation techniques in at least a large N approach [51]. Alternatively, one can study the Bethe-Salpeter equations for a bound state, something which has received widespread attention in the meson/baryon sector [52,53,54]. Nevertheless, similar approaches in the glueball sector are still hard to handle [13].
this expression can be easily expanded around p 2 = 0, where we defined the IR moments Notice that (19) defines a formal power series, so it does not need to converge.
By passing to we arrive at a system with finite boundaries.
Before going any further, there is a technical issue to deal with. We notice that eq. (11), or equivalently eq. (18), is a divergent integral, simply visible on dimensional ground through power counting. This is a typical feature of quantum field theory and the standard way to obtain a finite spectral integral is to subtract the first few orders of its Taylor expansion [55]. If we need r subtractions at p 2 = T , with T the subtraction scale, the subtracted spectral representation reads: We will thus consider the moment problem associated to the moments ν n , obtainable viâ Analogously as before, we introduce with The final question is, given a set of moments ν 0 , . . . ν 2N −1 , how to construct the underlying spectral density, with hopefully good convergence properties in terms of the numbers of input moments? For this so-called reduced Hausdorff problem, a nice answer has been provided in terms of Padé approximants of order [N , This is nothing else than a rational function approximation to the original function f (z ) since P i and Q i are polynomials of the designated order. The associated solution for the spectral density reads i.e. the spectral function is approximated by a (finite) series of δ-functions. This sounds well-suited for particle physics where peaks in the spectral functions correspond to particles 9 . The mass estimates, the t i 's, precisely correspond to the poles of the Padé approximant. We will not need the B N i 's, but they are pivotal to answer another important question: when does a set of moments correspond to a positive spectral function, i.e. a probability distribution? The answer is when the residues, i.e. the B N i 's, are positive numbers. We refer 9. δ-peaks correspond to stable particles, while finite peaks with a certain width to unstable particles. An inherent shortcoming of Padé approximation is thus that unstable particles are also replaced by δ-peaks. to [56,57,58,59] for details and proofs, including the interesting properties of the Padé approximants. We limit ourselves to mention here the link with orthogonal polynomials and the location of the poles w.r.t. the cut structure of f (z ): the Q N form a set of polynomials orthogonal over [0, s 0 ] with weight σ(s ). Consequently, their zeroes z * will all be real, different and lying in the interval ]0, s 0 [. Essentially, the poles of the rational approximant have replaced the branch cut of the original function. Now it is clear that any "reasonable approximative calculation scheme" (being the moment problem) capable of generating poles in the subtracted result sub (p 2 ), shall also give poles for (p 2 ), since the difference between both expression is just a polynomial in p 2 with infinite coefficients. Thus, we could equally well search for the poles of f (z ) in terms of z (≡ −1/p 2 ), which we know to be located in the interval ]0, s 0 [, viz. for p 2 ∈] − ∞, −1/s 0 [. An expansion at small p 2 corresponds to a Laurent expansion in 1 z near z ∼ ∞. By power counting, we see that f (z ) ∼ 1 z for z ∼ ∞, consistent with a [N , N − 1] Padé approximant. At lowest order, the Padé approximant will in general look like which in return will lead to a mass estimate m g l u e b a l l = ν 0 Returning to the question of how we are actually approximating the original (one-loop) correlation function: we start from a (necessarily subtracted) one-loop approximation to the composite operator's propagator using a nonperturbative input gluon propagator. Since in any glueball channel, it is, to our knowledge, quasi impossible to resum classes of diagrams using approximations as in NJL-models for quark bound states, it looks unfeasible to explicitly construct a pole in this correlator by resummation techniques 10 . Since we anyhow have at our disposal a first order approximation for the bound state correlation function, the proposed Padé approximation amounts to replace this first order function with a rational one that (i) does have a pole, and (ii) approximates the original function well. From the latter perspective, we point out that for N → ∞, the corresponding Padé approximant will converge to the original (one-loop) function f (z ) everywhere except on the branch cut of f (z ) where the poles of the approximant will pile up, with the smallest pole moving closer to the branch point of f (z ) [58]. Though, since the function f (z ) that we are Padé approximating is only a lowest order truncation of the full correlator, it would seem obsolete to consider a high order Padé approximant. To merely illustrate how "well" a lowest order Padé approximant performs, we have displayed in Figure 1 the function f (z ) defined in eq. (24) together with its approximation (28) for the SU (3), 0 ++ case. We recall that the approximation is in terms of large z and we see that the original function f (z ) -not displaying a pole-is quite well replaced by a rational function -with pole. The function f (z ) does not exist for 0 ≤ z ≤ s 0 as this corresponds to the original branch cut in momentum squared space. With the hindsight that Padé approximation does solve the Hausdorff moment problem, a lowest order Padé approximation does correspond to keeping only the first two   moments of the moment problem. One might wonder why it was a priori necessary to go through all the troubles of deriving the (one-loop) spectral representation of the correlation function, why not rather compute for Euclidean momenta p 2 ≥ 0 the (subtracted) one loop Feynman integral and directly 11 expand this up to order p 2 ? This recipe would also immediately be applicable to whatever gluon propagator one would like to use as input. Though, we believe this is the wrong way to proceed, since in that case also unphysical pieces of the, in some approximation computed, correlation function will sneak into the Padé approximation, e.g. from using a gluon propagator with cuts and/or poles in the complex plane as with the RGZ one or as with the fitting form proposed in [45]. With the current setup, we are assured that only the physical piece (positive spectral density, branch cut along the negative real axis) enter the rational approximation scheme.

(Pseudo)scalar and (pseudo)tensor glueball operator in d = 4
Having armed ourselves with the necessary technology, we are ready to present some results. First, we consider the first 2 moments corresponding to the spectral densities (12)- (15). Figure 2(a) shows the mass estimate in terms of the variable subtraction scale T for the gauge group SU (3), based on the RGZ fit (2),(3). A simple power counting argument reveals that the (pseudo)scalar spectral representation is well-defined (finite) after 3 subtractions, while the (pseudo)tensor sector needs 7 subtractions. As tried in [16,19], we could search for an optimal T in terms of which there is a minimal dependence of physical variables on the a priori free T , following the spirit of [60]. Unfortunately, there is no optimal T to be found. In fact, for all allowed values of T > −τ 0 there exists a mass estimate. Therefore, we suggest to sample over all possible T 's to get an average mass estimate. A logical choice seems to call for a Gaussian distribution centered at zero momentum subtraction (the latter being a rather common choice), whilst allowing for a variable width ω > 0. We will then fix ω using the principle of minimal sensitivity (PMS). Thus, It is clear from Figure 2(b) that there is a minimal ω-dependence in the scalar and pseudoscalar case. Albeit less clear from that same Figure 2(b), the (pseudo)tensor case displays an inflection point. After numerically 11. One could even think of first expanding and then performing the loop integral to further simply the computations at hand.    computing the optimal values, ω PMS (0 ++ , 0 −+ , 2 ++ , 2 −+ ) ≈ (1.32, 1.45, 2.82, 3.32) GeV 4 , the optimal mass averaged estimates are shown in Table 1, along some lattice and a Hamiltonian quasi-particle model values, taken from the quoted papers. Where necessary, we converted to GeV units by using the typical value σ = 0.44 GeV, notice that this is essentially the value for the string tension discussed in [62,63]. We simply picked the central value of the lattice papers to get an idea of their mass estimates. In order to get a rudimentary error estimate on the reported mean value m , we computed the standard deviation σ associated to the distribution (30) for ω = ω PMS , finding σ(0 ++ , 0 −+ , 2 ++ , 2 −+ ) ≈ (0.26, 0.27, 0.67, 0.79) GeV.

Scalar and tensor glueball operator in d = 3
As a final application, we consider 2 glueball states in d = 3, based on eq. (2) and eqns. (16)- (17). In d = 3, we need one less subtraction compared to d = 4, thus r = 2, resp. r = 6 for the scalar, resp. tensor glueball. The quoted analytical work [65] relies on the field correlator method, while [66] on the gauge invariant Karabali & Nair variables [67,68]. The corresponding numbers can be found in Table 3. Unfortunately, it is clear from Figures 4(a), 4(b) we can neither extract an optimal width ω nor an optimal subtraction scale T in the d = 3 case. At best, we can put an estimate 2-3.5 GeV with the 2 ++ probably heavier given its larger mass over most of the shown ω-(or T -) interval.

Discussion
We have discussed, with a recently introduced IR moment technique and ensuing rational (Padé) approximation, mass estimates for a set of glueball states in gluodynamics. The cases SU (3) and SU (2) have been reported in d = 4, while in d = 3 the case of SU (2) has been considered.
Our calculations have relied on the tree level confining gluon RGZ propagators, which encode nonperturbative information on the Gribov horizon and on dimension two condensates. Even if we have been working in a lowest order approximation, the fact that in most cases the Padé estimates for the masses reside in the same ballpark as other analytical and/or lattice estimates -which themselves sometimes show a widespread range of numbers-gives credit to our infrared moments method in general and to the expectation that important nonperturbative information is already present in the gluon propagator. Future work should be directed towards adding higher loop contributions, e.g. in the Refined Gribov-Zwanziger framework, to the glueball correlators and see how they influence the mass estimates, although it seems not easy to obtain the higher order corrections to the spectral densities in closed form. In such instance, the Padé method might also be of help since the first few moments of a specific correlator could be computed directly from the momentum space expression of the correlator and, as mentioned in the main text below eq. (27), the rational approximants can be used to decide whether these moments can belong to a positive spectral density or not [56]. It would also be instructive to investigate how sensitive the Padé results are to different numerical/functional prescriptions for gluon propagators, as obtained, for instance, in [46,45,70,71,72,73,74].