Spectrum of scalar and pseudoscalar glueballs from functional methods.

We provide results for the spectrum of scalar and pseudoscalar glueballs in pure Yang-Mills theory using a parameter-free fully self-contained truncation of Dyson-Schwinger and Bethe-Salpeter equations. The only input, the scale, is fixed by comparison with lattice calculations. We obtain ground state masses of 1.9 GeV and 2.6 GeV for the scalar and pseudoscalar glueballs, respectively, and 2.6 GeV and 3.9 GeV for the corresponding first excited states. This is in very good quantitative agreement with available lattice results. Furthermore, we predict masses for the second excited states at 3.7 GeV and 4.3 GeV . The quality of the results hinges crucially on the self-consistency of the employed input. The masses are independent of a specific choice for the infrared behavior of the ghost propagator providing further evidence that this only reflects a nonperturbative gauge completion.

Introduction

Glueballs, i.e. hadrons that consist of gluons only, are extremely fascinating objects to study. They arise due to the non-Abelian nature of Yang–Mills theory which allows for the formation of gauge invariant states of gluons that interact strongly amongst each other. The properties of glueballs have been studied in many models since their prediction in the 1970s [1, 2]. Today, many glueball masses in pure Yang–Mills theory are known rather accurately owing to high statistics quenched lattice calculations [3-6]. Unquenched lattice calculations of glueball masses are still on the exploratory level with considerable uncertainties due to severe problems with the signal to noise ratio, see, e.g. [7] and references therein. Alternative theoretical frameworks, such as Hamiltonian many body methods [8, 9] or chiral Lagrangians [10, 11], have shed some light on potential mass patterns and identifications of experimental states dominated by their glueball content, see [12] for a comprehensive review. However, it seems fair to state that our detailed understanding of glueball formation from the underlying dynamics of Yang–Mills theory is still far from complete. In this work, we provide an additional, complementary perspective from functional methods.

While the calculation of mesons and baryons from functional bound state equations is a very active field, see, e.g. [13, 14] and references therein, studies of exotic states in this framework are less abundant. This is particularly true for glueballs due to the inherent complexity of gauge-fixed Yang–Mills theories. Some results have been reported in [15-18], but they remain on an exploratory level due to the ansaetze used for the input. An alternative approach to extract glueball masses from Landau gauge correlation functions was followed in [19, 20].

Fig. 2

Interaction kernels from the three-loop 3PI effective action. All propagators are dressed; black disks represent dressed vertices. In our calculation, we include the diagrams inside the red rectangles. Details are discussed in the text

Functional methods allow direct studies of the internal structures of bound states as determined by the dynamics of QCD. As an example, let us mention tetraquarks where the role of different two-body clusterings can be studied, see Ref. [21] for a summary. Another promising aspect of functional methods is that they provide direct access to analytic properties. The necessary numerical techniques are not fully developed yet, but they have already been applied to some interesting questions like dynamical resonances or spectral functions, see, for instance, [22-30]. As a third example where functional methods can provide useful insights into QCD we mention the study of its phases and the transitions between them, see, e.g., [31, 32]. In the present work we add another application of functional equations to the list. It is set apart from the previous examples by delivering quantitative results without the need for any model parameters which typically substitute for missing information in the considered truncated systems of equations. This is to our knowledge the first such calculation. If such setups can be devised and solved also for other quantities of interest, this would constitute another major step forward.

What enabled us to perform a parameter-free calculation is the solid understanding of the properties of Yang–Mills correlation functions obtained in the last decade and a refinement of the methods to compute them from functional equations, e.g. [33-43]. In particular, we employ recent results from a fully self-contained truncation of Dyson–Schwinger equations (DSEs) for propagators and vertices that leads to good agreement with corresponding gauge fixed lattice results without any tuning [43]. As we will see in the following, the glueball spectrum extracted from the corresponding set of bound state Bethe–Salpeter equations (BSEs) agrees quantitatively with corresponding lattice results.

The remainder of the article is organized as follows. In Sect. 2 we introduce the bound state equations we solve and the employed input. The extraction of the spectrum is explained in Sect. 3 and the results are presented and discussed in Sect. 4. We summarize and provide an outlook in Sect. 5. The employed methods are illustrated in the appendix for a meson system.

Framework and input

We consider Landau gauge Yang–Mills theory with gauge fixing via the Faddeev–Popov procedure. The corresponding bound state equation for a two-gluon glueball is depicted in Fig. 1. Since for some quantum numbers the two-gluon state also couples to a ghost–anti-ghost state, we also need to consider the corresponding BSE. For simplicity, we refer to them as glueball-part and ghostball-part. The forms of the interaction kernels can be derived from an effective action and its truncation. We use the kernels obtained from a 3PI effective action truncated to three loops [44, 45]. This choice is motivated by excellent results that have been obtained in this truncation scheme for the coupled set of DSEs for the dressed gluon propagator, the dressed ghost propagator and the dressed ghost-gluon as well as three-gluon vertices [43] discussed below. A truncation on a similar level has also been used in the quark sector in Ref. [46]. The corresponding kernels, shown in Fig. 2, are derived by performing appropriate functional derivatives of the action [16, 47–49].

Fig. 1

The coupled set of BSEs for a glueball made from two gluons and a pair of Faddeev–Popov (anti-)ghosts. Wiggly lines denote dressed gluon propagators, dashed lines denote dressed ghost propagators. The gray boxes represent interaction kernels given in Fig. 2. The Bethe–Salpeter amplitudes of the glueball are denoted by gray disks

There are some noteworthy and important differences between a setup based on the 1PI effective action and our 3PI setup. The 3PI effective action depends on dressed propagators and three-point functions which can be calculated from their equations of motion corresponding to stationarity conditions of the action. Typically, one employs a (dressed) loop expansion for the 3PI effective action with a finite number of terms [44, 45]. The truncated DSE system derived therefrom is then a finite and closed system of equations for propagators and three-point functions with all higher-order vertices bare. Thus, it contains non-perturbative information for the three-point functions in a significantly different fashion than a truncation derived from the 1PI effective action. In the latter, all propagators and vertices are treated on the same footing, and for any truncation of the 1PI effective action, there are infinitely many equations to be solved. Thus, one normally performs the truncation at the level of the system of equations (derived from the 1PI effective action) instead of the action itself. If the 3PI effective action is truncated to three loops, the equations of motions for the three-point functions have a similar structure as the equations from the 1PI effective action without two-loop terms [44]. However, there are no bare three-point functions in the former case. This seems like a small difference, but it turned out that the equations of motion of the three-loop expanded 3PI effective action, which are one-loop, outperform the results from equations of 1PI equations truncated to one loop. An explicit comparison for the three-gluon vertex showed that including the two-loop diagrams for the latter leads to good agreement. However, such two-loop calculations are much harder to realize [43].

The similar structure of the equations of motion from 1PI and 3PI effective actions persists also for the kernels of the bound state equations. A typical kernel for a 1PI calculation contains only one-particle exchanges with one bare vertex, but some studies beyond that exist as well, e.g., [50, 51]. The lowest order of 3PI kernels contain only one-particle exchanges but in contrast to the 1PI case all vertices are dressed. The next order also contains two-particle exchanges. For now, these are not taken into account due to their technical complexity.

Of course, a priori there is no known way to guarantee that higher order corrections are small in a dressed loop expansion. The results of this work, in particular the excellent agreement of our glueball masses with corresponding lattice results, provides some circumstantial evidence that this may be the case. Further evidence has been collected in [35, 36, 38, 40], where selected extensions of the present truncation scheme have already been calculated and found to be small indeed. Nevertheless, more systematic tests in dedicated calculations need to be done in the future to further elaborate on this issue.

Having specified the underlying equations, we can now turn to solving them. A solution to the coupled set of BSEs for glueballs can be found by treating them as an eigenvalue equation for the matrix of interaction kernels with the eigenvector combining the glueball-part and the ghostball-part of the Bethe–Salpeter amplitude. This leads towhere P and p are the total and relative momenta. For a solution, the eigenvalue must equal 1. The mass is determined from the corresponding value of the total momentum, .

This structure is very general and valid for all allowed quantum numbers. Here we treat the simplest glueballs, the scalar and pseudoscalar ones. To this end, we need to specify the forms for the amplitudes. The glueball-part has two open Lorentz indices and is transverse. For the scalar glueball () one can find two independent tensors with these properties given by [15]:with . For the pseudoscalar glueball (), only one tensor exists which we choose asThe hat indicates normalization and the superscript T that the vector is made transverse with respect to P. For the ghostball-part, which is a scalar in Lorentz space, the amplitude is simply given byThere is no corresponding amplitude with negative parity. This simplifies the BSE for the pseudoscalar glueball where the ghostball-part of the amplitude drops out.

The input required to solve the BSEs are the dressed gluon and ghost propagators, and , respectively, given byas well as the dressed three-gluon and ghost-gluon vertices. For these we use numerical results from a DSE system also derived from the three-loop 3PI effective action. A graphical representation of this truncation together with a thorough discussion of all technical details and merits can be found in [43]. Here, we only wish to mention that the scheme is self-contained, i.e., it can be solved without any ad-hoc ansaetze and parameters. Thus, either correlation functions are taken into account and solved for self-consistently, or they are consistently neglected.

A truncation independent property of the Yang–Mills correlation functions in the continuum is that they appear as a one-parameter family of so-called decoupling/scaling solutions with the scaling limit as an end-point [56-59]. For the ghost and gluon propagators, these solutions agree at large momenta but start to differ around the scale of and below. The appearance of this one-parameter family has been discussed to be caused by incomplete nonperturbative gauge fixing in the Landau gauge/Faddeev–Popov setup due the Gribov problem [58, 60]. Indeed, on the lattice it is well studied that different samplings of the Gribov region influence the low momentum behavior of correlation functions, see, e.g. [60-62] and references therein. If this picture is correct, all solutions of the family should lead to the same results for physical quantities. The glueballs studied in this work provide an ideal testing ground for this hypothesis.

We thus used a number of different solutions within this one-parameter family including the scaling solution taken from Ref. [43] for our calculation of glueball masses. The corresponding gluon propagators and dressing functions are shown in Fig. 3. Since the precise correspondence of continuum and lattice gauge fixing methods is an open issue [53], we show a selection of functional results and lattice results for one particular prescription to deal with Gribov copies. All employed gluon propagators feature a maximum at , but they are so shallow for two cases that they are hardly visible in the plots. The ghost dressing function is shown in Fig. 4. For the plots, the propagators were renormalized such that they agree with the lattice data of Ref. [52] at some high momentum. This rescaling is unrelated to the renormalization employed in the actual calculation, where different renormalization conditions were used. Especially for the vertices the consistency of the renormalization procedure is an important issue as discussed in [43]. Here, we can directly use these renormalized results.

Fig. 3

Gluon dressing function (left) and gluon propagator (right) in comparison to lattice data [52]. For the sake of comparison, the functional results were renormalized to agree with the lattice results at . Different lines correspond to different decoupling/scaling solutions as explained in the text

Fig. 4

Left: Ghost dressing function in comparison to lattice data [52]. For the sake of comparison, the functional results were renormalized to agree with the lattice results at . Right: Ghost-gluon vertex dressing function at the symmetric point in comparison to SU(2) lattice data [53]. Different lines correspond to different decoupling/scaling solutions as explained in the text

The ghost-gluon vertex dressing function is shown in Fig. 4. The plot only shows the symmetric point at which all momentum squares are equal, but its full kinematic dependence was taken into account in the calculation. The ghost-gluon vertex is finite in the Landau gauge [63] and thus no renormalization is necessary. The three-gluon vertex is shown in Fig. 5 where it is compared to lattice data from Refs. [54, 55]; see Refs. [64, 65] for similar results. Also here the full kinematic dependence was taken into account, but its angular dependence is very weak. For the sake of comparison, all data were renormalized to 1 at .

Fig. 5

Solutions for the three-gluon vertex dressing function at the symmetric point in comparison to lattice data [54, 55]). For the sake of comparison, all data were renormalized to 1 at . Different lines correspond to different decoupling/scaling solutions as explained in the text

With the functional renormalization group, a truncation very similar to the one used here leads to equally good results for the propagators and vertices [37] which further supports the trustworthiness of the present truncation. While improvements of this truncation were partially already tested and found to be small, e.g., the impact of neglected contributions like four-point functions [36, 38, 40] and a full basis for the three-gluon vertex [35], more tests should be performed in the future.

For quark anti-quark bound states, the consistency of the truncations for the DSE and BSE system is known to be crucial. This is true in particular for the chiral properties, where the axial-vector Ward–Takahashi identity relates the self-energies and integration kernels, see, e.g. [14]. In the present case of pure Yang–Mills glueballs, similar relations do not exist to our knowledge, but it is certainly important to consider the impact of consistency between the truncations of the DSEs and the BSEs.

Extraction of the spectrum

Bethe–Salpeter equations describe bound states with time-like total momenta. As a consequence, propagators and vertices inside Bethe–Salpeter equations are tested at complex squared momenta. While calculations in the complex plane are in principle possible with functional methods, see, e.g. [22, 25, 27, 30, 66–69], such calculations are realized only with comparatively simple truncations. More advanced truncation schemes, like the one we use in this work, have been solved for space-like Euclidean momenta only. Analytic continuations of correlation functions have been explored [70-73], but definite conclusions are still lacking.

Thus, we follow here an alternative approach and instead analytically continue the eigenvalues calculated for Euclidean momenta into the time-like momentum domain . For every bound state we have calculated the eigenvalue of the BSE for 100 space-like and real momenta and extrapolate to time-like momenta using Schlessinger’s method based on continued fractions [74, 75]. The extrapolation error is dominated by (small) numerical inaccuracies in the input data and the limitations of the chosen extrapolation function itself. We partially quantify the extrapolation error by a bootstrap-like procedure: We take only a random subset of 80 points, which is enough to get a good extrapolation function, and calculate the mass. Excluding exceptional extrapolations that do not give a bound state, we repeat that process 100 times and average the results. The one standard deviation errors indicated in our results for the bound state masses stem from this procedure. Of course, this does not include additional errors due to the truncation and other sources.

We tested the extrapolation procedure in a case where we are able to compare directly with a solution. To this end, we used a quark–anti-quark BSE in rainbow-ladder approximation that is reviewed in Ref. [14]. The comparison between the extrapolated and the calculated eigenvalues is described in Appendix A. Within errors both results agree very well.

One complication which needs to be taken into account is that not every eigenvalue curve in the space-like domain can be extrapolated to a physical state with a positive real mass. We encountered (and discarded) some complex eigenvalue curves and one example for a (spurious) tachyonic state in case of the pseudoscalar glueball. Furthermore, eigenvalue curves may cross. Hence, the hierarchy of eigenvalues at space-like momenta does not need to correspond to the hierarchy of the solutions.

The extrapolation is illustrated in Fig. 6 where the solutions for off-shell values of are shown. We chose deliberately to print the eigenvalues over and not M, as the extrapolation is done in . For this plot, Eq. (1) was solved for various values of . As explained above, the errors are calculated by averaging over 100 random subsets of the input points. The Schlessinger extrapolation sometimes contains artifacts, e.g., very narrow poles. These can affect the average locally which can lead to the small roughness observed in the errors. It might be possible to refine the averaged extrapolations by applying further statistical methods or other optimization techniques as, for example, in Ref. [72]. For now, the obtained precision is sufficient.

Fig. 6

Eigenvalues for off-shell values of for scalar (left) and pseudoscalar (right) glueballs. Shown are the eigenvalues for the ground state and the first two excited states

Finally, the scalar glueball BSE leads to some additional insight with regard to the ultraviolet behavior of the Bethe-Salpeter amplitudes. The integral in the BSE is convergent only if the amplitude falls off polynomially at large momenta. This is a standard behavior of all Bethe–Salpeter amplitudes in the quark sector [14] and also expected for the glueball one. However, we do find eigenvalue curves related to states with amplitudes that run logarithmically and therefore do not satisfy this constraint.1 As a consequence, these eigenvalue curves are cut-off dependent. In the gluon DSE, a similar problem is associated with terms that break gauge covariance and lead to quadratic divergences. A range of methods have been applied to remove these terms, see e.g. [76] and references therein. In the BSE, the simplest way to deal with this problem is to identify the solutions with the wrong asymptotics as artifacts and discard them. We adopted this strategy in this work for both scalar and pseudoscalar glueballs. It remains to be seen whether the problem disappears entirely in even more advanced truncation schemes.

Results and discussion

We show our results for the ground and excited states of scalar and pseudoscalar glueball masses in Table 1 and Fig. 7 together with results from lattice calculations. After the first version of this article appeared on arxiv.org, new lattice data became available [6] which we also include now in the comparison. For the sake of comparability, we adjusted all scales to the same value of the Sommer parameter . Originally, we set the scale of our results from the gluon propagator of Ref. [52] where was used, while for the lattice results of Refs. [4, 5] it was . In Ref. [6], was set to to express the obtained results in physical units. Here, we adopt this value for the comparison of all results and rescaled the masses accordingly. For our results this means that we multiply the masses by the factor 1.059. In addition, we show the results in terms of the corresponding scalar ground state masses.

Table 1

Ground and excited state masses M of scalar and pseudoscalar glueballs. Compared are lattice results from [4–6] with the results of this work. For [4, 5], the errors are the combined errors from statistics and the use of an anisotropic lattices. For [6], the error is statistical only. In our results, the error comes from the extrapolation method. All results use the same value for , see text for details. The related error is not included in the table. Masses with are conjectured to be the second excited states

State	[4]		[5]		[6]		This work
State	M[MeV]	M/M_0^++	M[MeV]	M/M_0^++	M[MeV]	M/M_0^++	M[MeV]	M/M_0^++
0^++	1760(50)	1(0.04)	1740(50)	1(0.04)	1651(23)	1(0.02)	1850(130)	1(0.1)
0^^∗++	2720(180)	1.54(0.11)	–	–	2840(40)	1.72(0.034)	2570(210)	1.39(0.15)
0^^∗∗++	–	–	–	–	3650(60)^†	2.21(0.05)^†	3720(160)	2.01(0.16)
					3580(150)^†	2.17(0.1)^†
0^-+	2640(40)	1.5(0.05)	2610(40)	1.5(0.05)	2600(40)	1.574(0.032)	2580(180)	1.39(0.14)
0^^∗-+	3710(60)	2.1(0.07)	–	–	3540(80)	2.14(0.06)	3870(120)	2.09(0.16)
0^^∗∗-+	–	–	–	–	4450(140)^†	2.7(0.09)^†	4340(200)	2.34(0.19)
					4540(120)^†	2.75(0.08)^†

Fig. 7

Results for scalar and pseudoscalar glueball ground states and excited states from lattice simulations [4, 6] and this work. In the left diagram we display the glueball masses on an absolute scale set by . In the right diagram we display the spectrum relative to the ground state

For the scalar glueball we find a ground state mass of . The extrapolation is very stable in this case. The first excited state is at , and we also find a candidate for the second excited state at . The pseudoscalar glueball ground state is at . The first two excited states are at and . The masses for the scalar and pseudoscalar second excited states were, to our knowledge, predicted here for the first time.

Ref. [6] contains lists of states for the different irreducible representations of the octahedral subgroup to which the rotational symmetry is broken on the lattice. The spin of a glueball state is identified by searching for nearly degenerate states in corresponding representations. The scalar and pseudoscalar glueballs only appear in the representation . Four states with the corresponding quantum numbers , are listed. The lower two are identified as the (pseudo)scalar glueball and its first excitation. The other two are very close to each other. Since their values are similar to the ones found here for the second excited states, we conjecture that one of each pair indeed is such a state.

The comparison with the lattice results of Refs. [4-6] shows that we do not only recover the same hierarchy but that our results also agree quantitatively very well with the lattice results even including second excited states.

The amplitudes for the ground states and the two excited states at the lowest calculated value for are shown in Fig. 8. Although they do not correspond to the physical amplitudes, for which , they provide a lot of insight. First of all, we see for the meson example described in detail in Appendix A that the amplitudes show only a small dependence on the momentum variable . In that case, we can even use the results at space-like momenta to extrapolate the results to time-like momenta. Second, off-shell amplitudes are required for calculations where the corresponding bound states are intermediate states over which one integrates.

Fig. 8

Leading amplitudes for the scalar (top) and pseudoscalar (bottom) ground states (left), first excited states (middle) and second excited states (right) for the lowest calculated

All amplitudes in the ground states show a distinct peak around . The first and second radial excitations contain one and two nodes, respectively. For the scalar glueball we observe an interesting interplay of the three different components given in Eqs. (2) and (4). Whereas the ghostball-part plays an important role in the ground state, it is negligible for the two excited states. Thus, neglecting the ghost contribution would distort the results quantitatively (but not qualitatively) as we tested explicitly. Furthermore, for ground and excited states different gluon amplitudes are largest.

We explicitly verified that the masses extracted for the glueballs are invariant (within error bars) of our choice of input within the one-parameter family shown in Figs. 3, 4 and 5. The endpoint of this family is the scaling solution for which we obtain also the same results. This supports the view that the family corresponds to different gauge choices within Landau gauge [60]. At the same time, it constitutes a nontrivial test of the reliability of the employed truncation, as any inconsistency could destroy the gauge independence easily. We explicitly verified that the spectrum gets distorted if inconsistent input (e.g. propagators from one member of the family and vertices from another) is used. This sensitivity on the employed input indicates that it would be challenging to replace part or all of the input by models.

In this context, it is also interesting to mention that in the process of devising the setup employed in this work, we also tried various models for propagators and/or vertices for testing purposes. However, not only were we never able to produce results somewhat close to the present ones, also qualitative features of the eigenvalue spectrum were different. As an example, let us mention Ref. [16], in which, using model vertices, a reasonable result for the scalar glueball was obtained, but the pseudoscalar one was above . We thus conclude that the quality of the input is very important for a good overall picture.

Summary and outlook

We calculated the masses of the three lowest states for scalar and pseudoscalar glueballs in pure Yang–Mills theory from propagators and vertices in the Landau gauge obtained from Dyson-Schwinger equations. We gave special emphasis to the consistency between the DSE and BSE setups. The truncation we employed is completely self-contained and does not depend on any parameter or ansatz. The scale is inherited from the comparison with lattice results for the gluon propagator. The gauge dependent input for the propagators and vertices translates into invariant results for the masses within a one-parameter family of nonperturbative completions of Landau gauge. Our results agree quantitatively very well with the lattice results of Refs. [4, 5] and we add two more states to the known spectrum of pure Yang–Mills theory. It is encouraging that these new states can be matched with subsequently published lattice results [6].

To our mind, the present results constitute a considerable step forward to describe QCD bound states in a continuum approach from first principles. There are many applications to which this setup can and should be extended. In particular, the quark sector should be included to get access to real-world glueballs. The calculation of higher spin glueballs or even other exotics like hybrids may be possible within the same framework as well. At the same time, these results motivate further studies of functional equations for complex momenta to eliminate the need of the extrapolation of the eigenvalue curves from space-like Euclidean data. A step in this direction was taken in Refs. [25, 30].