Higher spin glueballs from functional methods

We calculate the glueball spectrum for spin up to J=4 and positive charge parity in pure Yang-Mills theory. We construct the full bases for J=0,1,2,3,4 and discuss the relation to gauge invariant operators. Using a fully self-contained truncation of Dyson-Schwinger equations as input, we obtain ground states and first and second excited states from extrapolations of the eigenvalue curves. Where available, we find good quantitative agreement with lattice results.


Introduction
Hadrons that strictly consist of only gluons may not exist in nature. Instead, one expects mixtures of these socalled glueballs with corresponding meson states with the same quantum numbers. Nevertheless, it is important to study the spectrum of glueballs in pure Yang-Mills theory, since it may very well be that some of these states obtain only small corrections from the matter sector of QCD. In this respect, it is very interesting to note that the scalar glueball predicted by lattice Yang-Mills theory long ago [1,2,3,4] and recently also within functional methods [5] seems to show up in radiative J/Ψ -decays [6] with almost unchanged mass. This merits further investigation in approaches that can deal with the (anti-)quark admixtures of full QCD such as unquenched lattice calculations [7], the functional a e-mail: markus.huber@physik.jlug.de b e-mail: christian.fischer@theo.physik.uni-giessen.de c e-mail: helios.sanchis-alepuz@silicon-austria.com approach [5,8,9,10,11,12,13], Hamiltonian many body methods [14,15] or chiral Lagrangians [16,17], see also [18,19,20,21,22] for reviews.
In a previous work [5], we provided first results for the masses of ground and excited glueball states in the scalar and pseudoscalar channels of pure Yang-Mills theory using a functional approach based on a fully selfcontained truncation of Dyson-Schwinger equations and a set of Bethe-Salpeter equations derived from a threeparticle-irreducible (3PI) effective action. In this work, we generalize our framework and present additional results for ground and excited states with quantum numbers J = 2, 3, 4 and positive as well as negative parity.
The obtained results are completely model independent in the sense that there are no free parameters which need to be tuned. For the input, the only truncation appears at the level of the 3PI effective action which is truncated to three loops and allows a self-consistent determination of all the required correlation functions. This feature of our calculation is different from other functional calculations which rely on phenomenological motivated model interactions. Conceptually, it is appealing since extensions like the inclusion of quarks are naturally contained in this framework. Technically, it comes at the price of more demanding calculations to obtain the input.
The manuscript is organized as follows. We first introduce the bound state equation in the next section, followed by a discussion of quantum numbers and how to construct bases for the Bethe-Salpeter amplitudes in Sec. 3. The input and the procedure to obtain the glueball spectrum are explained in Sec. 4 and the results are presented in Sec. 5. We conclude with a summary in Sec. 6.

Glueball bound state equations
The properties of particles, elementary or composite, can be extracted from correlation functions. We describe here how this is done for glueballs and how the Bethe-Salpeter formalism is formally related to the calculation of glueballs from gauge invariant operators on the lattice. For pure (anti-)quark bound states, this is discussed, e.g., in [23]. For glueballs, we have to pay particular attention to the fact that individual twoor three-gluon states themselves cannot be gauge invariant. This is obvious when considering the operators used in lattice calculations to calculate glueball masses. On the other hand, the functional formalism provides a means to extract the gauge invariant mass of bound states from gauge variant n-point functions. This is particularly convenient for glueballs, where we can then get such information from the two-body bound state equations alone. The properties of a particle can be extracted from an appropriate two-point function. As a prime example for glueballs let us consider the composite operator O 1 (x) = F µν (x)F µν (x). It has the quantum numbers J PC = 0 ++ [24]. In momentum space, the two-point function D(P ), has a pole at the mass of the 0 ++ glueball. This is, for example, exploited in lattice calculations where the mass of a particle can be extracted from the exponential decay of such a two-point correlator at large times: Since the operator O 1 is gauge invariant, also the twopoint function and the masses are gauge invariant. Expanding the operator in gluon fields, one can write down a functional equation [25]. In the present example, this leads to an expression with up to six loops where all quantities are fully dressed. The derivation is fully worked out for a similar example in Ref. [26] and can also be done algorithmically with a computer [27]. O 1 can be split into contributions with two, three, and four gluon fields, and one can group their contributions to D(P ) accordingly into terms with four, five, . . . , eight gluon fields. A crucial point in determining the mass of a state is that its pole already shows up in one or more of the individual contributions, which therefore may be studied separately.
The one-loop contribution to D(P ) stems entirely from the two-gluon contribution to O 1 , i.e. the correlation function of four gluons G µνρσ :  Fig. 1, couples a contribution containing two gluons (which we call glueball-part) and one containing a ghost and an antighost (ghostball-part). The ghost fields appear in the gauge fixing procedure; here we work in Landau gauge. The corresponding kernels in the BSE have been derived from the 3PI effective action expanded to three loops [28,29,30,31,32] and are shown in Fig. 2. The derivation is detailed in Ref. [5]. The necessary input, i.e. the two and three-body correlators appearing in the kernels, is also obtained from this effective action [33]. In our main analysis, we only take into account the diagrams leading to one-loop expressions in the BSE. These are enclosed in the red boxes. The inclusion of two-loop diagrams has to be deferred to future work, as they are computationally much more expensive. For now, we take the quantitative agreement of our oneloop results with lattice results as indication that these diagrams are dominant. Note also that decays are not possible in this setup as it would require additional diagrams, see, e.g., [34].

Quantum numbers and tensor bases
In this section, we explain which quantum numbers we calculate and derive the corresponding tensor bases. = + = + 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. 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.

Landau-Yang theorem
Landau and Yang determined the allowed values for spin and parity quantum numbers of a system consisting of two photons [35,36]. They concluded that J = 1 in general and P = −1 for odd spin are forbidden. This is often directly transferred to glueballs 'consisting of two gluons'. We emphasize that one cannot construct a gauge invariant operator from two gluon fields and in this sense no pure two-gluon glueballs exist. Here we consider the two-gluon part of the full gauge-invariant equation which can contain a pole by itself. Below we discuss why the Landau-Yang theorem does not apply here and does not remove the pole from the two-gluon part.
The Landau-Yang theorem in the context of non-Abelian gauge theories was already considered in Refs. [37,38,39,40]. The decisive difference of QCD to QED is the color degree of freedom which leads to more allowed quantum numbers for color antisymmetric states. However, such states are not relevant for the case of two gluons in a color singlet state.
Of direct relevance to our framework is one particular assumption that enters in the derivation of the Landau-Yang theorem, namely that the photons/gluons are onshell. The derivation proceeds by considering the parity even and odd cases separately. In the former case, the state J P = 1 + is removed from the possible states by the on-shell condition. For the parity odd case, the states J P = (2k + 1) − , with k a nonnegative integer, are removed also by the on-shell condition and the fact that the spatial dimension is three [35,39]. Since the gluons in the glueball are not on-shell, the on-shell condition is violated and the J = 1 states are not removed.
Neither are the ones with odd J and negative parity. Thus, the Landau-Yang theorem does not forbid glueballs with these quantum numbers in our formalism.
Whether corresponding states actually appear as solutions of the two-body BSE is entirely a dynamical question. We will see below, that in the current truncation scheme this is not the case.

Charge parity
As described above, we employ a two-body bound state equation to calculate glueball masses. This gives access to states with charge conjugation parity C = +1 as we explain below. Since the gluon field itself carries a color charge, it transforms nontrivially under charge conjugation. The SU (3) color matrices λ a get transposed under a charge conjugation. Based on their symmetry properties, the gluon field transforms as [41] A a where η(a) = +1 a = 1, 3, 4, 6, 8 −1 a = 2, 5, 7.
A color neutral operator with two gluon fields thus transforms as This corresponds to positive charge parity. Negative charge parity, on the other hand, can be realized, for example, with three gluon fields: The last equality follows from the fact that the only nonvanishing elements of the symmetric structure constant d abc have zero or two indices equal to 2, 5 or 7.
For the antisymmetric structure constant, only the elements with one or three indices equal to 2, 5 or 7 are nonzero. Consequently, this would lead to positive parity.

Spin and parity
For the tensor bases we need to find all tensors compatible with a given spin J and parity P. In addition, the transversality of the gluon propagator in Landau gauge needs to be taken into account. We describe how to construct such tensor bases and work them out explicitly for J = 0, 1, 2, 3, 4. For glueballs, we only need to consider integer spin states. Spin J can then be described by a tensor T of rank J. It must be symmetric, traceless and transverse with respect to the total momentum P [42]: The second and third conditions remove lower spin states from the tensor and we end up with 2J + 1 independent elements [42]. These properties can be enforced on any tensor using appropriate spin projection operators P J , see, e.g., [43]. Let us consider as an example the case of spin two for which the projection operator reads P J=2 ρσρ σ (P ) = 1 2 (P ρρ (P )P σσ (P ) + P ρσ (P )P σρ (P )) where is the transverse projector. Applying P J=2 on any tensor T ρ σ , the result obeys the conditions (8): It is symmetric in the two indices ρ and σ, the trace is zero as ensured by the last term and it is transverse with respect to the total momentum P . For the spin-2 example we can now easily construct a basis for the ghostball-part. We need tensors with two Lorentz indices. Using the metric tensor, the total momentum P and the relative momentum p we obtain five different structures. The transversality condition (8c) leaves us with only p µ p ν and g µν . The tracelessness condition (8b) further eliminates the latter. Thus, the only admissible tensor for the spin-2 ghostball-part is For other spins, the ghostball-part tensors are constructed in the same way by projecting J relative momenta with a spin projector. For the glueball-part, the procedure is similar, but now two additional Lorentz indices from the gluon legs enter the game. For reference, we call them 'gluon leg indices' in contrast to the 'spin indices' ρ and σ. It is thus advantageous to filter the starting set of tensors appropriately, because for higher spin J, the number of possible tensors with J + 2 indices increases quickly. In the case of spin two, we can reduce the number of tensors constructed from the metric g and two momenta from 43 to 10. First of all, many tensors will not survive the spin projection due to the transversality or tracelessness conditions. This can easily be taken into account by discarding all tensors with a spin index attached to a total momentum or metric tensors in the spin indices. Furthermore, the symmetry condition entails that upon spin projection many tensors become linearly dependent. To avoid that, we only take one representative for each case. Following these selection rules, we determine for spin two the following initial set of tensors: {g µρ g νσ , g µν p ρ p σ , g µρ p ν p σ , g νρ p µ p σ , p µ p ν p ρ p σ ; g µρ P ν p σ , g νρ P µ p σ , p µ P ν p ρ p σ , P µ p ν p ρ p σ , P µ P ν p ρ p σ }.
µ and ν are the indices of the gluon legs and ρ and σ are spin indices. We call this set a 'pre-basis'. As a next step, we consider the transversality of the gluon propagator in the Landau gauge. It entails that the gluon legs of the glueball-part are transversely projected on the right-hand side of the BSE. To be consistent with the left-hand side, the BSE amplitude should not contain nontransverse parts as otherwise the equation is no longer a proper eigenvalue equation. Consequently, we need to modify the pre-basis such that it is invariant under P(p 1 ) µµ P(p 2 ) νν , where p 1 and p 2 are the momenta of the gluon legs, see Fig. 3 for our naming conventions. The gluon momenta are related to the total and relative momenta by Below we will use all four momenta P, p, p 1 , p 2 for convenience, although of course only two of them are linearly independent. A simple consequence of the transverse projection from the gluon propagators is that the momenta P and p attached to a gluon leg index are no longer independent and the pre-basis is reduced further. In the example of spin two, one can choose the first five tensors in Eq. (12). These still need to be 'transversalized'. We illustrate this step explicitly in Sec. 3.4 for spin zero. Roughly speaking, it amounts to replacing momenta with a gluon leg index by a transversalized momentum. In our numerical treatment of the glueball BSE we find that the use of a transverse base is crucial to obtain correct results. At this point we have what we call a 'transverse prebasis'. To obtain the final bases, the spin projection operators P J are applied to the elements of the prebases. This is discussed explicitly in Secs. 3.4-3.8. For negative parity states we additionally need the Levi-Civita symbol ε which requires a special discussion. Due to its antisymmetric nature, it can appear only in a limited number of ways and we list the possible variants. First, it can have either of the gluon indices µ or ν. Second, it can have only one index from the spin part because of the symmetric nature of the spin projector. Third, it can have indices contracted with the relative or total momentum. This leads to the following possibilities: As stated above, any permutations of the spin indices are irrelevant because of the symmetry of the spin projector. We now list the possible tensors for positive and negative parity. In the former case, we take all tensors which can be constructed from the momenta and the metric g. We do not need to consider an even number of Levi-Civita symbols ε, as such an expression is linearly dependent on the previously mentioned terms. This can be directly seen by writing the product of two Levi-Civita symbols in terms of metric tensors.
The following list of possible tensors for positive parity is exhaustive. Further required indices are filled by ap-propriate expressions as discussed below. The symbol '•' denotes a Lorentz index from the spin part.
-Two metric tensors: There is only one tensor with g µ• g ν• . -One metric tensor: There are three possibilities, namely g µν , g µ• or g ν• .
-No metric tensor: There is exactly one tensor constructed from momenta only.
At this point we can already count the number of tensors for each J. To this end, we observe that spin indices need to come with the relative momentum p due to the transversality of the spin projector. Open gluon leg indices, on the other hand, can come with two different momenta which, however, are linearly dependent when transversely projected. Thus, every tensor will be represented only once. We obtain two tensors for J = 0 (of the types g µν and p µ p ν ), four for J = 1 (the type g µ• g ν• does not exist) and five for higher spin (all possibilities realized).
The tensors ε µραβ p α P β P ν , ε νραβ p α P β P µ are not included, because they are linearly dependent after transverse projection of the gluon legs. Counting the possible numbers of tensors taking into account the transversality of the gluon propagators leads to one for J = 0 (of the type ε µναβ p α P β ), five for J = 1 (of the types ε µναβ p α P β , ε µνρα p α , ε µνρα P α , ε µραβ p α P β p ν , ε νραβ p α P β p µ ) and seven for higher spin. However, these numbers reduce as some tensors are linearly dependent. This can be shown by an explicit calculation. The final numbers are then three for J = 1 and four for higher spin. We note that there is no ghostball-part for negative parity, because in that case we only have spin indices at our disposal of which we cannot put more than one into the Levi-Civita symbol. We will now list the pre-bases for the glueball-parts for different spins constructed from the arguments above. We also symmetrize the expressions with respect to exchanging the two legs. For all basis tensors, the spin can be inferred from the number of indices. Tensors with/without tilde refer to positive/negative parity.

J = 0 glueballs
The scalar and pseudoscalar glueballs naturally have the simplest basis. We are going to illustrate the step from the initial basis to the pre-basis, the transversalization, with the scalar one. For positive parity, one can write down five tensors following the arguments outlined above: Taking into account the transversality of the gluon propagator, the last four tensors are equivalent, since P(p 1 ) µµ P µ ∝ P(p 1 ) µµ p µ and similarly for the other leg with index ν. There are, thus, only two tensors in the transverse part. Naively, one could take the tensors from the first line, but they are not invariant under a transverse projection of the gluon legs. One could achieve this by simply projecting transversely. Instead, we are going to use the auxiliary tensor where i, j = 1, 2. For i = 1, it is transverse in p 1,µ and similar for j = 2 in p 2,ν . Two transverse tensors constructed from t ij µν are then An appropriate normalization was chosen to make the tensors dimensionless. This basis from Ref. [10] was employed in Ref. [5]. 1 We also tried the variant based on t 11 µα t 22 αν , but it led to numerical instabilities. The ghostball-part for the scalar glueball is simply a scalar. For negative parity, only one tensor exists which can be chosen as In Ref. [5] we used The hat indicates normalization and the superscript T that the vector is made transverse with respect to P . However, the two expressions are equivalent except for the normalization, as also the first tensor is transverse in the gluon legs. This can be seen from the antisymmetry property of the Levi-Civita symbol due to which it can only hold two linearly independent momenta. From the transverse projectors of the gluon propagators thus only the metric part survives. 1 We corrected here a typo for the momenta with index ν.

J = 1 glueballs
For spin one, we construct a basis from the set of all tensors with three Lorentz indices. As we saw for the scalar glueball, it is sufficient to consider only the relative momentum p for the gluon leg indices due to the linear dependence introduced by the transverse projection. For positive parity, we obtain then the tensors g µν p ρ , g µρ p ν , g νρ p µ , p µ p ν p ρ . We recombine them such that they are symmetric under exchange of the gluons and invariant under transverse projection of the gluon legs. Furthermore, we normalize them conveniently. First, however, it is useful to introduce the following auxiliary momenta: Our choice for the transverse pre-basis is To single out the contributions relevant for spin one, we apply the spin-1 projector, which is identical to the transverse projector with respect to the total momentum: For negative parity, the following tensors can be constructed: After spin-1 projection, only three of these tensors are linearly independent. The first three tensors are a possible choice for the pre-basis.

J = 2 glueballs
The initial tensors for spin two were discussed in Sec. 3.3. Following the transversalization procedure outlined above, we obtain the following transverse pre-basis: where An appropriate normalization of all tensors was chosen for convenience. Applying the spin-2 projector (9) yields the basis. The possible tensors for negative parity are Of these seven tensors, four are linearly independent after transverse projection. We choose the tensors τ 1 , τ 2 , τ 3 , and τ 6 .

Composite operators for glueballs
The bases we constructed can be related to the physical picture of glueballs described in form of gauge invariant operators as given, e.g., in Ref. [24]. Let us consider the following examples: Evidently, these operators contain two-, three and fourgluon components. As discussed in Sec. 2, poles in the correlation functions of any of these components correspond to glueball masses. Since we consider a twogluon bound state equation, we expect that the bases for the Bethe-Salpeter amplitudes have a correspondence in the two-gluon component of such gauge invariant operators. This can be checked by applying two derivatives with respect to gluon fields and setting the fields to zero. We exemplify this with the operator T , from which we obtain 10 terms: g µν g ρσ →τ 1 µν g ρσ (J = 0) g µρ g νσ , g µσ g νρ →τ 1 µνρσ (J = 2) g µν p 1ρ p 2σ , g µν p 2ρ p 1σ →τ 2 µνρσ (J = 2) g µρ p 1ν p 2σ , g µσ p 1ν p 2ρ →τ 3,4 µνρσ (J = 2) g νρ p 2µ p 1σ , g νσ p 2µ p 1ρ →τ 3,4 µνρσ (J = 2) g ρσ p 2µ p 1ν →τ 2 µν g ρσ (J = 0). (34a) The first and last, which vanish upon spin-2 projection, are identified with a basis for the scalar glueball multiplied by g ρσ , see Sec. 3.4. Making these tensors transverse, relates them to the basis in Eq. (18). Some of the other tensors are identical upon projection with the spin-2 projector. For instance, this applies to the second and third tensors. Upon closer inspection, we find that tensors two and three correspond (after transversalization) to τ 1 µνρσ , tensors four and five to τ 2 µνρσ , and tensors six to nine to τ 3 µνρσ and τ 4 µνρσ of Eq. (25). The third column in (34) indicates these correspondences. The basis element τ 5 µνρσ does not appear here. It is interesting to note that we find the amplitude of this tensor to be subleading.

Solving the BSEs: input and extrapolation
We follow the same setup as in [5] which we shortly summarize here. To solve the BSE, we need the ghost and gluon propagators and the ghost-gluon and three-gluon vertices. We take them from a parameter-free calculation [33]. The corresponding quantities are compared to lattice results in Figs. 4 and 5. For a meaningful comparison, the functional results were renormalized to the values indicated in the captions. For the calculations, the original data was used which was renormalized as explained in Ref. [33]. They also agree very well with results from the functional renormalization group [51]. We emphasize that a family of solutions can be obtained which all fulfill the standard Landau gauge fixing condition [52,53,54,55,56,57,58,59]. These solutions may correspond to different nonperturbative gauge completions of the perturbative Landau gauge [55]. In Ref. [5] we tested the dependence of our results for J = 0 on the different solutions. We found that for all solutions, including the one with a diverging ghost dressing function, that the obtained masses agree within the errors induced by the extrapolation. Hence, we continue here with only one solution, the one shown in Figs. 4 and 5. It is worth mentioning that these results are independent of the number of colors as the system scales trivially with g 2 N c . For Yang-Mills theory, nontrivial terms only enter at four loops [60]. The truncation in [5] does not capture these terms. As the employed BSE also scales trivially in g 2 N c , our results are hence independent of N c . Since the two-and three-point functions necessary as input for the BSEs are only available for Euclidean momenta, we cannot solve the BSE directly at the timelike pole-locations, −M 2 = P 2 < 0. To access time-like momenta, we instead extrapolate the eigenvalue curve from the space-like side using Schlessinger's method based on continued fractions [61,62]. We tested the reliability of the eigenvalue extrapolation by calculating a meson with a setup that allows calculations for timelike momenta [5]. In that case, the method was very accurate for masses up to roughly 2 GeV. Beyond that, the errors increase. For the present case, where no direct estimation of the extrapolation error is possible, we estimate it indirectly by taking 100 samples of 80 points out of 100 calculated points for the extrapolation. From this we calculate the error indicated in the results section. The error estimated from this method establishes only a lower bound and results for large masses should be interpreted with care. The identification of states is also somewhat involved for several reasons. First, it is not necessary that the hierarchy of eigenvalues at space-like momenta agrees Gluon and ghost dressing functions Z(p 2 ) and G(p 2 ), respectively, (left) and gluon propagator D(p 2 ) (right) in comparison to lattice data [45]. For the sake of comparison, the functional results were renormalized (in the plots only) to agree with the lattice results at 8 (ghost) and 6 GeV (gluon).  [46]. Right: Three-gluon vertex dressing function at the symmetric point in comparison to lattice data [47,48]), see Refs. [49,50] for similar results. For the sake of comparison, the three-gluon vertex data was renormalized (in the plot only) to 1 at 5 GeV.
with the hierarchy in the physical region, since eigenvalue curves may cross at time-like P 2 . It is then not clear, whether the extrapolation method is still reliable. Second, we also encountered cases where it was not possible to identify a unique eigenvalue curve with a state, because there are several solutions leading to similar masses. We indicate all such cases in the discussion in the next section. For higher spin states the numerical effort to determine excited states increases drastically due to the larger number of Lorentz indices. As a consequence we had to reduce the numerical accuracy in order to comply with the available CPU time. We checked the effect of this on the scalar and pseudoscalar glueballs and found that the ground states are unaffected but excited states may disappear with decreased accuracy. This explains why we do not always find two excited states for J > 0, as discussed below. For spin three and four and posi-tive parity, we also did not take into account the diagrams containing ghosts to reduce the computational cost. However, we confirmed for one eigenvalue that the amplitude from the ghostball-part is subleading compared to the ones from the glueball-part.

Results
Our results for glueball states with J = 0, . . . , 4 and positive and negative parity are summarized in Fig. 6 and Tab. 1. For the comparison with lattice results, we rescale all results to the same value r 0 = 0.472(5) fm of the Sommer scale used in [4], see Ref. [5] for details. The results for J = 0 are taken from Ref. [5] and are unambiguous. We found clear signals for the ground states and two excited states which compare well with lattice results. Note that the second excited states on [2] [  [2,3,4] with the results of this work and [5]. For [2,3], the errors are the combined errors from statistics and the use of anisotropic lattices. For [4], the error is statistical only. In our results, the error comes from the extrapolation method and should be considered a lower bound on errors. All results use the same value for r 0 = 1/(418(5) MeV). The related error is not included in the table. Masses with † are conjectured to be the second excited states. Masses with * come with some uncertainty in their identification in the lattice case or in the trustworthiness of the extrapolated value in the BSE case.
the lattice were not uniquely identified in Ref. [4]. In fact, for P = +1 and for P = −1 two states were found in each case with very similar masses, see Tab. 1. The good agreement with our results indicate that one of them is most likely indeed a second excited J = 0 state.
The identification of states for J = 2 is more challenging than for spin zero. For 2 ++ , we observed a degeneracy in the ground state and first excited masses, viz., we found two eigenvalue curves for each that yield, within extrapolation errors, the same masses. The curves themselves, however, are not degenerate. Inspecting the Bethe-Salpeter amplitudes from each of these two solutions, we identified a crucial difference: one set does not show the power law fall-off in the large momentum region as expected for a normalizable Bethe-Salpeter amplitude. We adopted this as a criterion to rule out solutions of this type, as discussed already in Ref. [5]. Within error bars, the two remaining solutions for the ground and first excited states are in very good agreement with corresponding lattice results. We also found a candidate for a second excited state at 3.8 GeV, but the error from the extrapolation is more than 1 GeV in that case so we refrain from including it in the table.
Our results for the masses in the 2 −+ channel again show an interesting pattern with seemingly two solution branches. Inspecting the amplitudes of one branch we find a clear pattern with respect to zero crossings: no crossings for the ground state and n zero-crossings for the n-th excited state. The amplitudes furthermore resemble qualitatively those of the pseudoscalar glueball. This pattern is not seen, however, for the second branch of solutions, which we therefore discard as artifact of the truncation. In the physical branch, we find an excited state in the mass region beyond 4 GeV. Since its extraction involves quite some extrapolation into the time-like region, its mass value comes with a large uncertainty.
For glueballs of spins three and four, we obtain masses between 3 and 4.5 GeV. Again, we assign some additional uncertainty to these values due to the extrapolation. Let us first discuss the positive parity states. For 3 ++ , we find two excited states. All three states are rather close to each other. It remains to be seen whether this still persists in an improved approach where we would be able to follow the eigenvalue curves into the time-like region. In principle, we cannot exclude that two of the curves merge at some point and disappear BSE lattice [Morningstar , Peardon , 1999 ] lattice [Athenodorou , Teper , 2020 ] Fig. 6 Results for glueball ground states and excited states for the indicated quantum numbers from lattice simulations [2,4] and functional equations. In the upper plot, we display the glueball masses on an absolute scale set by r 0 = 1/(418(5) MeV).
In the lower plot, we display the spectrum relative to the ground state. Masses with † are conjectured to be the second excited states. Masses with * come with some uncertainty in their identification in the lattice case or in the trustworthiness of the extrapolated value in the BSE case.
(i.e. the eigenvalues become complex). Such a behaviour has already been seen in approaches with simpler truncations, where time-like momenta are accessible, see e.g. [23,63] for an overview. In our case, this would mean that not all three states would survive. In contrast, our ground state for 4 ++ can be identified exceptionally cleanly and is in the same range as the corresponding lattice result. The extrapolation is extremely stable in that case what leads to the unusually small error. The reason for that is currently unknown but we also observed it for some other eigenvalues. As far as lattice results are concerned, we also have to keep in mind that the identification of the 3 ++ and 4 ++ states is not as unique as for the lighter states [4]. The negative parity states 3 −+ and 4 −+ have not been identified on the lattice and we regard our results as predictions, albeit with large uncertainties, especially since the spin four state is lighter than the spin three state. Finally, we also solved the BSE for spin one. As discussed in section 3.1, there is no principal reason why such a state should not appear in a two-body bound state equation. However, dimensional counting of the corresponding operators containing these quantum numbers suggests that the masses should be large [2,64]. Indeed, in the mass region considered in this work we were not able to find solutions fulfilling the criteria for the amplitudes discussed above. For 1 ++ , this agrees with the results from lattice Yang-Mills theory, while for 1 −+ a mass of about 4 GeV was found [4].

Summary and discussion
We calculated the spectrum of glueballs with C = 1 and J = 0, 1, 2, 3, 4 from a two-body bound state equation. For J = 1, we did not find solutions. For the other cases, we were able to find the ground states and up to two excited states. The results for J = 0 have been discussed already previously in Ref. [5]. A crucial factor to arrive at these results was the use of high-quality input for the gluon and ghost correlation functions from a self-contained calculation [33]. As the only parameter of the input is the scale, this is a calculation from first principles. Our results agree qualitatively and quantitatively with corresponding lattice results [2,3,4], especially when comparing the ratios to the lightest state. This is very encouraging. The largest deviations we observe for the 2 −+ and 3 ++ groundstates. The present setup can be improved in several ways. An obvious one concerns the inclusion of the neglected twoloop diagrams in the BSE to establish its self-consistency. This is computationally very costly but in principle possible. However, due to the good agreement with lattice results so far, we do not expect large corrections at least for the lighter states. More importantly, we rely on an extrapolation of the eigenvalue curves from space-like to time-like momenta. A calculation directly at P 2 = −M 2 could shed light on some of the uncertainties we encountered in the identification of physical states and remove the extrapolation error. Unfortunately, for such a calculation we require two-and threepoint functions evaluated in the time-like (and complex) momentum regime, which is currently not available in sufficient quality. Some corresponding calculations exist but are less advanced in their truncations [65,66,67]. Of course, the most obvious and interesting extension of the present framework concerns the inclusion of the matter sector of QCD. This is work in progress.