Yang-Mills glueball masses from spectral reconstruction

We compute masses of the two lightest glueballs from spectral reconstructions of timelike interaction channels of the four-gluon vertex in Landau gauge Yang-Mills theory. The Euclidean spacelike dressings of the vertex are calculated with the functional renormalisation group. For the spectral reconstruction of these Euclidean data, we employ Gaussian process regression. The glueball resonances can be identified straightforwardly and we obtain $m_{sc} = 1870(75)~$ MeV as well as $m_{ps} = 2700(120)~$ MeV, in accordance with functional bound state and lattice calculations.


I. INTRODUCTION
The hadronic spectrum of Yang-Mills theory and QCD includes purely gluonic bound state contributions, the glueballs.The experimental verification of their existence is an important test of QCD; however, it is not yet conclusive [1][2][3][4] as these states are difficult to access due to their large overlap with other hadronic resonances.Different possible experimental candidates have been proposed, including the various f 0 states, some of which are expected to appear in decay channels of J/ψ [5,6].The overlap with other states also complicates their theoretical determination when considering QCD; for corresponding lattice calculations see [7][8][9].In Yang-Mills theory, the situation is much simpler and the first few lightest states are well known; for lattice results see e.g.[8,[10][11][12][13][14].For computations with functional approaches-in particular with a combination of Dyson-Schwinger equations (DSE) and Bethe-Salpeter equations (BSE)-see e.g.[15][16][17][18][19][20].
In this work, we put forward a self-consistent functional ansatz for computing masses of bound states by exploiting their overlap with resonant interaction channels of gauge-fixed correlation functions.The approach is then used to determine the masses of the scalar (J P C = 0 ++ ) and pseudo-scalar (J P C = 0 −+ ) Yang-Mills glueballs, utilising the fact that these states have overlap with channels of the four-gluon vertex that carry the respective symmetries, where they appear as peaks of the corresponding spectral functions.We use Gaussian process regression (GPR) to compute these spectral functions by reconstructing Euclidean correlators obtained within the functional renormalisation group (fRG) framework in [21].The inversion of the spectral representation is an ill-conditioned problem; see e.g.[22][23][24].The applicability of GPR to such linear inverse problems was discussed in [25] and the approach has since been employed to compute ghost and gluon spectral functions from 2+1 flavour lattice QCD results [26].* j.turnwald@theorie.ikp.physik.tu-darmstadt.de This paper is organised as follows.In Section II, we introduce the spectral representation of Euclidean dressing functions and discuss the projections onto the fourpoint vertices in Yang-Mills theory.The reconstruction approach using GPR is described in Section III.In Section IV, the resulting spectral functions are presented and we report the masses of the scalar and pseudo-scalar glueballs.We conclude in Section V.

II. SPECTRAL REPRESENTATIONS OF YANG-MILLS CORRELATION FUNCTIONS
Non-perturbative calculations of correlation functions in Yang-Mills theory are generally only possible in Euclidean space-time, either on the lattice or with functional approaches.While the latter framework in principle also allows direct access to real-time propertiesalbeit with a qualitatively increased effort-real-time lattice calculations are often faced with intractable signalto-noise problems.Accordingly, computing timelike observables such as transport coefficients, pole masses, and decay rates typically requires the reconstruction of timelike correlation functions from their spacelike Euclidean counterparts via the associated spectral representations.

A. Spectral representations
Correlation functions of physical states, and in particular the two-point functions, admit a spectral representation.For the propagator-the inverse 1PI two-point function-this is the Källén-Lehmann (KL) representation, Here, G(p) denotes the Euclidean propagator and ρ(ω) the spectral function, which is obtained by The spectral functions of asymptotic states are positive semi-definite and admit the interpretation of a probability density.In gauge theories, however, the situation becomes more complicated.To begin with, even the existence of a KL representation is not settled for ghost and gluon propagators, and (1) may feature additional structures in the complex momentum plane; for a detailed discussion see [21].Moreover, in the Landau gauge, the gluon and ghost spectral functions exhibit negative infrared (IR) and ultraviolet (UV) tails; see [27].These properties can be inferred from the respective IR and UV asymptotic behaviour of the Euclidean correlation functions [21,[27][28][29], and these relations also hold true for the present analysis involving four-gluon vertices.Note also that while gauge-fixed correlation functions may not permit a KL representation, the scattering matrix elements are directly constructed in terms of these correlators and obey (1).Hence, features which are in direct correspondence to observables-such as bound states-can still be extracted from such gauge-fixed correlation functions.

B. Four-gluon correlation function
In the present work, we consider single interaction channels that have overlap with the bound states of interest.The spectral representations of these channels follow directly from the structure of the full, analytically continued correlation functions; see e.g.[30].Up to minor modifications, they are given by (1): for the relevant scalar dressings of the four-gluon vertex, we use [31] ( The constant part λ A 4 ,0 accounts for the classical contribution.The Euclidean dressings of the interaction channels are computed with the fRG; for a recent review see [32].The diagrammatic representation of the associated equation is shown in Figure 1; more details on the fRG approach and the specific computation for the vertex are deferred to Appendix A. We remark in this context that correlation functions in Landau gauge Yang-Mills theory computed within sophisticated truncations to the fRG pass all available lattice benchmark tests; see [21,33].This concerns in par-ticular the ghost and gluon propagators, whereas lattice results for vertices still exhibit large uncertainties.Nevertheless, since state-of-the-art functional results for correlation functions fully agree with lattice calculations within statistical errors, any reconstruction based on the former approaches is consistent with the latter.
In order to access the masses of the scalar J P C = 0 ++ and pseudo-scalar J P C = 0 −+ glueball, we have to determine tensor structures and momentum channels that overlap with these states.In general, it is desirable that the chosen channels have overlap only with the states of interest, as any reconstruction method faces increasing problems with multi-peak structures due to the exponential suppression of heavier states in the Euclidean data.Accordingly, their resolution requires an exponentially increasing accuracy, contributing to the ill-conditioned nature of the reconstruction problem.
For the scalar glueball, the above requirement is particularly simple to satisfy, since it is the lightest excitation and the classical tensor structure suffices, i.e.
Correspondingly, we use for the pseudo-scalar glueball, which has the correct transformation properties (see e.g.[15,19]), and does not overlap with the scalar glueball.In ( 5), ε µνρσ denotes the fully antisymmetric tensor and the momenta are chosen to be orthogonal, p 1 • p 2 = 0. Finally, we have to specify the momentum channels: we restrict ourselves to a single exchange momentum and the external (incoming and outgoing) momenta are chosen to have the same magnitude, p 2 ≡ p 2 1 = p 2 2 .This leaves us with two invariants: p 2 and x = (p 1 • p 2 )/p 2 .For the scalar glueball, the momenta are chosen to be parallel, i.e. x = 1; for the pseudo-scalar one, they are chosen to be orthogonal, i.e. x = 0. Further details on the projection operators are given in Appendix A 2.

III. GAUSSIAN PROCESS REGRESSION WITH INDIRECT OBSERVATIONS
GPR is widely employed as a non-parametric interpolation method for noisy observations.In essence, GPs can be used to define probability distributions over families of functions that fit a given set of data without explicitly assuming a functional basis.For an in-depth introduction to GP theory and applications, see e.g.[34].
Recently, GPR has also been applied to the probabilistic inversion of the KL spectral representation [26] as well as the extraction of parton distribution functions [35][36][37].The present work follows the same line of reasoning: making use of the fact that GPs are closed under linear transformations, it is possible to infer data from indirect observations that are related to the quantity of interest by a linear forward process [25].In particular, we may obtain predictions for the spectral function from measurements of the associated correlator without inverting the KL transformation directly.
To this end, we start by defining a GP prior distribution over spectral functions that encodes our knowledge and assumptions about ρ(ω) before making any observations, Evaluating this GP for any set of points ω i results in a multivariate normal distribution with mean µ(ω i ) and covariance k(ω i , ω j ).Furthermore, the discrete propagator data G(p i ) are then also normally distributed, with mean and covariance obtained by applying the linear forward process (the KL integral) to µ and k, i.e.
As GPs can be specified completely by their second-order statistics, µ is usually set to zero for simplicity since any information contained therein may be fully absorbed into k.However, a non-zero prior mean may still be useful in practice, in which case it can simply be subtracted from the data beforehand.Using bold symbols for vectors of discrete data, e.g.p for a set of momenta p i , the joint distribution of spectral function values ρ at any point ω and a set of correlator data G(p) can be expressed as with Since the joint distribution is normal, the posterior distribution of the spectral function conditioned on observations of the correlator can be derived in closed form, This is a standard result in multivariate statistics and is essentially equivalent to GPR with direct observations, only with additional insertions of the linear transformation that one seeks to invert.The GP posterior (10) encodes our knowledge of the spectral function given the correlator data and directly accounts for some additive Gaussian noise with variance σ 2 n in the observations.For computational applications of GPs, the covariance is usually parameterised by a kernel function k(ω, ω ′ ), as already implied by the notation used above.Since this fully specifies the GP as mentioned previously, choosing the right type of kernel is a pivotal part of finding a good model for a given set of data.A natural choice in many applications are so-called universal kernels that can describe any continuous function [38] and hence provide the required flexibility when little is known about other properties of the desired solution a priori.The radial basis function (RBF) kernel, is a popular choice due to its universality and every function in the associated prior being infinitely differentiable, and is also employed in the present work.The parameters ℓ and σ control the length scale and overall magnitude of the correlation between data and are subject to optimisation; see Appendix B.
Predictions obtained with GPR can also be understood within the well-known Backus-Gilbert framework [39], one of the most popular approaches to spectral reconstruction in the lattice community.In fact, both methods produce numerically equivalent estimates under certain conditions [25], despite following different philoso-phies.Nevertheless, the GPR picture is much more flexible, since essentially any available prior information can be systematically incorporated into the regression by extending the covariance matrix in (8), following the same reasoning as in the construction of the joint distribution of G and ρ.Such prior information may simply consist of known values of the spectral function at certain points, in which case they are treated as direct observations.More generally, it can be any indirect data related to ρ through a linear operator-such as a derivative [40]-and even inequality constraints such as bounds and monotonicity conditions [41].
In summary, GPR is a powerful approach to tackle ill-conditioned linear inverse problems probabilistically, which makes it an attractive candidate algorithm for spectral reconstruction in quantum field theory.

IV. RESULTS
We calculate the vertex dressings with the fRG as outlined in Section II B; for details on the truncation and computation, see Appendix A. The resulting Euclidean dressing functions for the scalar and pseudo-scalar projections are shown in Figures 2a and 3a, respectively.
In the channels considered here, the ghost loops drop out; see Appendix A 2. Hence, these channels are free of the IR divergences that are in general present in the fourgluon vertex and we can utilise the constraint ρ A 4 (0) = 0 in the GP reconstruction.Furthermore, an additional bias is introduced in order to suppress unphysical oscillations at the tails of the spectral function.Similar to the procedure applied in [26], this is achieved by rescaling the frequency with a soft step function, where the parameter ℓ 0 controls the steepness and ω 0 the position of the midpoint.This rescaling can be understood as the introduction of a frequency-dependent length scale in the RBF kernel, with smaller values around ω 0 and larger values at the tails of the spectral function.We note that the resonances of interest are already observed without introducing this additional bias.However, the peaks are enhanced by this procedure while the reconstruction of the correlator remains in good agreement with the input data.While this parameterisation suppresses additional structures (such as excited glueball states at higher energies; see e.g.[19]), even without the rescaling (12) no additional features beyond the dominant peak corresponding to the bound state are observed, apart from the usual oscillatory behaviour at the tail of the spectral function.This implies that higher excited states exhibit at most sub-leading contributions to these vertex projections.Resolving these structures therefore requires either more sophisticated projections of the tensor structures or a significantly higher precision in the calculation of the vertex itself.The parameters of the RBF kernel and frequency rescaling are optimised by minimising an objective function, conventionally taken to be the negative loglikelihood (NLL).Unsurprisingly, the NLL shows a flat direction where some parameters are unconstrained; see Figure 4.This can be interpreted as a manifestation of the ill-conditioned nature of the inverse problem, and may be treated by imposing a hyperprior.We observe that changing the parameters in this direction has negligible impact on the resonant peak position; see Figures 5  and 6.Hence, the seemingly heuristic use of a generic hyperprior is well justified in this context as it does not introduce a bias for the quantity of interest.Details about this procedure as well as the optimised parameter values are provided in Appendix B. The intrinsic error estimate of the GP posterior is fixed to σ n = 10 −2 , corresponding to an upper bound on the uncertainty of the fRG calculation.σ n is not optimised as this diminishes the significance of the likelihood for the other parameters [42].
The results of [10,12] are rescaled to match [13,14] with r0 = 1/418(5) MeV.The errors of [10,12] are a combination of statistical as well as systematic uncertainties stemming from the lattice anisotropy and the scale r0.The errors for [13] are statistical only.For [14], the quoted values are the statistical as well as systematic uncertainties for the continuum extrapolation, respectively.For [19], the error comes from the extrapolation method.
menta, the result deviates more strongly, in particular for the scalar glueball.This is due to the additional bias introduced to the kernel that specifically suppresses any dynamics in the UV regime.The glueball masses are extracted from the dominant peak positions of the spectral functions.We obtain ωsc = 0.93 GeV for the scalar and ωps = 1.35 GeV for the pseudo-scalar channel.Since we work within the schannel approximation and have two incoming momenta each with the magnitude p, the peak position corresponds to the half of the glueball mass, i.e. m sc/ps = 2 ωsc/ps .Hence, we obtain the masses m sc = 1870(75) MeV for the scalar and m ps = 2700(120) MeV for the pseudo-scalar glueball.The reported errors are a combination of the standard deviations computed from the GP posterior and an additional 3% error from the scale setting procedure of the input data.A more in-depth discussion of the systematic error of the reconstruction can be found in Appendix B. We compare our results with masses obtained from independent lattice and DSE/BSE studies of the glueball spectrum in Table I and find them to be in reasonable agreement, in particular for the pseudo-scalar channel where they match well within the provided uncertainties.

V. CONCLUSION
We put forward a self-consistent approach for the extraction of bound state information from gauge-fixed correlation functions.Key to this framework is the spectral reconstruction of interaction channels in Euclidean spacetime that have overlap with the corresponding gauge-invariant bound state.The method is applied to lowlying glueball states in Yang-Mills theory, extracted from the dressing functions of the Euclidean four-gluon vertex.With appropriate projection operators of the four-gluon vertex, we obtain access to the masses of the scalar and pseudo-scalar glueballs.
The Euclidean dressings are obtained with the functional renormalisation group, also utilising earlier results for correlation functions from [21].The respective spectral functions are then computed via Gaussian process regression and their resonance peaks are identified with the glueball masses: for the scalar and pseudo-scalar glueballs, we arrive at 1870(75) MeV and 2700(120) MeV, respectively.The results agree well with independent studies of the glueball spectrum, lending further credibility to our proposed method of computing bound state properties from vertex dressing functions via spectral reconstruction.The present approach can also be directly applied to higher glueball states in Yang-Mills theory, as well as glueball and other hadronic states in QCD.

fRG equation for the four-gluon vertex
The master equation of the fRG is the flow equation of the scale-dependent 1PI effective action.It is obtained by introducing an IR cutoff with a cutoff scale k via a momentum-dependent mass function R k (p 2 ) that is added to the inverse propagator.The respective flow equation is derived via taking a derivative of the gener-ating functions w.r.t. the cutoff scale k, where t is the RG-time, and the trace in (A1) sums over species of fields, space-time (momentum), Lorentz indices and group indices.The regulator functions carry the classical dispersion of the ghost and gluon fields as well as a dimensionless shape function.The present results are computed with the usual exponential shape function, and an additional wave function renormalisation Z A,k or Z c,k ; for more details see [21].For a recent review of the fRG see e.g.[32] and references therein.Our general setup in Landau gauge Yang-Mills theory follows [21,33].The flow of the four-point vertex is obtained by taking the fourth derivative of (A1) w.r.t. the gluon field.In this work, we are only interested in certain channels of the four-gluon vertex.Hence, we do not solve the full system self-consistently, but take all other correlation functions such as the gluon propagators from [21] as input.
The fRG equation for the four-gluon vertex solved in the present work is depicted in Figure 1.This flow is integrated on the solution of the correlation functions obtained in [21].There, different IR closures of correlation functions in the Landau gauge have been computed, and the present work utilises the scaling solution.The independence of this choice has recently been shown in [19], where both solutions-decoupling and scaling-were considered in the context of glueballs.The approximation used in [21] only includes the primitively divergent (classical) tensor structures, which leads to semi-quantitative results.Further details can be found in [21].
We use the k-dependent dressing functions from [21] as input.Their parameterisations are given by Γ (2),ab ,abc A 4 ,µνρσ (p 1 , p 2 , p 3 ) = λ A 4 (p) f abn f cdn δ µρ δ νσ + perm., (A3) where we approximate the full momentum dependence of the vertices with the symmetric point configuration p, see e.g.[33], defined by with n = 3, 4. FIG. 4. Grid scans of the NLL (B1) of the reconstructions for both channels.Note that the optimisations are performed subsequently, starting with the RBF and followed by the bias parameters.The red lines indicate the trajectories in parameter space used for comparing the variance in the spectral functions; see Figures 5 and 6.The red cross indicates the NLL optimised parameters; see Table II.

Glueball projection operators
The full projection operator for obtaining the scalar glueball mass is simply a contraction with the transvere part of the classical tensor structure, given by (A5) Indices are suppressed for simplicity and the external momenta are already matched to the momentum parameter-isation of the four-gluon vertex, (A6) The classical four-gluon tensor structure τ A 4 ,cl is given in (4) and the transverse projection operator is The pseudo-scalar projection operator is defined by the tensor structure (5) and given by where the in-and outgoing external momenta are chosen to be orthogonal,

(A9)
We note that the projection onto the ghost loop part of the flow analytically vanishes for both projections on the momentum configurations under consideration.This was observed for a similar momentum configuration in [43].
nullify the bias, such as ℓ 0 becoming large.Hence, the parameters are first optimised only considering the bare RBF kernel in order to obtain baseline values.Subsequently, the bias is introduced and its parameters are optimised given the RBF kernel calculated beforehand.This way, the position and size of the dynamical part of the spectral function are also subject to optimisation.
The parameters are optimised by performing a highresolution grid scan; see Figure 4. Their optimal values are provided in Table II.In the direction of the magnitude parameter σ, the NLL does not change significantly for larger values.Similarly, in [26] this parameter was observed to exhibit an open direction towards infinity and a hyperprior had to be introduced, which is a manifestation of the ill-conditioning of the inversion.The dependence of the spectral function on σ is plotted in Figures 5  and 6, showing that while it does impact the magnitude of the dominating peak, its position and other general features of the spectral function remain stable.Accordingly, the overall magnitude of the computed spectral functions should be taken with a grain of salt, but predictions of other features such as the peak position, width, and overall shape are robust as the NLL diverges quickly when considering non-optimal parameters.Scanning the spectral functions in the plane of the bias parameters on the other hand reveals a more drastic change in the peak position.However, these parameters are restricted to a much smaller region by the likelihood and the stability of the peak position is retained.For a more quantitative statement about the systematic error of the reconstruction, an empirical comparison of different bias parameterisations is required.This can potentially be achieved by mapping out the posterior probability landscape with Monte Carlo methods.

FIG. 1 .
FIG. 1. fRG equation for the s-channel four-gluon vertex dressing.Wiggly orange lines correspond to fully dressed gluon propagators; black dots indicate fully dressed vertices.Permutations include the various possible configurations of external legs as well as permutations of the regulator insertion (indicated by a crossed circle).

FIG. 2 .
FIG. 2. (a)Euclidean dressing of the four-gluon vertex λ A 4 ,sc with the projection to obtain the scalar glueball mass, see Section II B, from the fRG (black crosses).This is compared to the reconstruction from the GP (green line).The corresponding spectral function ρ A 4 ,sc over frequency ω obtained with GPR is shown in (b).The light green band represents the 1σ region.
FIG. 3. (a)Euclidean dressing of the four-gluon vertex λ A 4 ,ps with the projection to obtain the pseudo-scalar glueball mass, see Section II B, from the fRG (black crosses).This is compared to the reconstruction from the GP (green line).The corresponding spectral function ρ A 4 ,ps over frequency ω obtained with GPR is shown in (b).The light green band represents the 1σ region.

( a )FIG. 5 .
FIG.5.Spectral function of the scalar channel.The bands show the variance with respect to the flat directions of the parameter space; see Figures4a and 4b.The peak position is observed to be robust, even under large parameter changes.The overall magnitude on the other hand shows considerable variation for both sets of parameters.