Dispersive analysis of glueball masses

We develop an inverse matrix method to solve for resonance masses from a dispersion relation obeyed by a correlation function. Given the operator product expansion (OPE) of a correlation function in the deep Euclidean region, we obtain the nonperturbative spectral density, which exhibits resonance structures naturally. The value of the gluon condensate in the OPE is fixed by producing the $\rho$ meson mass in the formalism, and then input into the dispersion relations for the scalar, pseudoscalar and tensor glueballs. It is shown that the low-energy limit of the correlation function for the scalar glueball, derived from the spectral density, discriminates the lattice estimate for the triple-gluon condensate from the single-instanton estimate. The spectral densities for the scalar and pseudoscalar glueballs reveal a double-peak structure: the peak located at lower mass implies that the $f_0(500)$ and $f_0(980)$ ($\eta$ ad $\eta'$) mesons contain small amount of gluonium components, and should be included into scalar (pseudoscalar) mixing frameworks. Another peak determines the scalar (pseudoscalar) glueball mass around 1.50 (1.75) GeV with a broad width about 200 MeV, suggesting that the $f_0(1370)$, $f_0(1500)$ and $f_0(1710)$ ($\eta(1760)$) mesons are the glue-rich states. We also predict the topological susceptability $\chi_t^{1/4}=75$-78 MeV, deduced from the correlation function for the pseudoscalar glueball at zero momentum. Our analysis gives no resonance solution for the tensor glueball, which may be attributed to the insufficient nonperturbative condensate information in the currently available OPE.


I. INTRODUCTION
We proposed to handle QCD sum rules [1] for nonpertrbative observables, like the ρ meson mass, as an inverse problem recently [2]. The spectral density on the hadron side of a dispersion relation, including both resonance and continuum contributions, is regarded as an unknown. The operator product expansion (OPE) on the quark side is calculated in the standard way. The spectral density (similar to a source distribution) is then determined from the OPE input (similar to a potential observed outside the distribution). This approach does not involve a continuum threshold, because the spectral density is supposed to be a smooth distribution, and may differ from the perturbative one (the quark-hadron duality is likely to be broken). It does not require a Borel transformation to suppress the continuum contribution, which will be solved from the inverse problem. The suppression on higher power corrections to the OPE can be achieved by considering the input in the deep Euclidean region. Once a dispersion relation is solved directly, the optional stability criterion [3,4] for sum rules is not necessary either. The long-existing concern on the rigorousness and predictive power of QCD sum rules [5,6] is then resolved. As an example, we demonstrated how to extract the masses and decay constants of the series of ρ resonances from the dispersion relation obeyed by a two-current correlator [2]. This new viewpoint, based only on the analyticity of physical observables, has been extended to the explanation of the D meson mixing parameters [7] and to the constraint on the hadronic vacuum polarization contribution to the muon anomalous magnetic moment [8].
The spectral density was parametrized as a sum of a pole contribution, which depends on the ρ meson mass and decay constant, and an arbitrary continuum contribution in [2]. The latter was expanded in terms of a set of orthogonal polynomials with unknown coefficients, which, together with the mass and decay constant, were deduced from the best fit of the hadron side to the OPE on the quark side. The continuum contribution turned out to be a ramp function, instead of a step function postulated in sum rules. A weakness of the above treatment is that the existence of the bound state, presumed from the beginning, is not a direct consequence of solving the inverse problem. In fact, the Bayesian approach [6] was applied to sum rules with a similar motivation, in which the form of a spectral density was not specified, but obtained using the maximum entropy method. It was found [9] that the resultant spectral density associated with the ρ meson exhibits a resonance peak at low energy, followed by a smooth continuum tail at high energy. Unfortunately, this approach is numerically intricate, and the peak location, i.e., the ρ meson mass, depends on some parameters introduced in the maximum entropy method. The complexity is expected, since it is notoriously difficult to solve an ill-posed inverse problem.
We will develop a simple inverse matrix method to solve for a spectral density from a dispersion relation without assuming the existence of a resonance: the whole spectral density is expanded in terms of the generalized Laguerre polynomials, whose coefficients are derived from the OPE input directly. The quark-hadron duality for a spectral density on the hadron side and the discretionary criteria for balancing perturbative and nonperturbative contributions [10,11] are not implemented. That is, unambiguous predictions for nonperturbative observables can be acquired in the above dispersive formalism. Besides, the precision of our predictions can be improved systematically by adding higher-order and higher-power corrections to the OPE input. We first test this proposal by solving for the ρ meson mass, and demonstrate that a resonance peak shows up in the spectral density naturally, a result the same as from the maximum entropy method [9]. It means that our formalism is equivalent to but technically simpler than the Bayesian approach [6]. By producing the ρ meson mass m ρ = 0.78 GeV [12] from the peak location, we fix the value of the gluon condensate α s G 2 = 0.08 GeV 4 , α s being the strong coupling constant, which is within the range accepted in the literature (see [13], for example). Our solution contains only a single peak corresponding to the ground-state ρ meson, manifesting the difficulty to explore properties of excited states with denser spectra, which is also encountered by the Bayesian approach [9].
The framework with the fixed gluon condensate is then applied to the analyses of the scalar, pseudoscalar and tensor glueball masses. Though the instanton background may contribute to the OPE [14], we will not consider its effect as in [13] to avoid the model dependence from, e.g., parametrizations for the instanton size distribution. We will show that the gluon condensates of various dimensions are enough for establishing the gluonium states, and confront our findings with those in [13]. For sum-rule investigations on scalar and pseudoscalar glueball properties in the instanton background, refer to [11,[15][16][17][18][19]. The quark-loop and quark-condensate corrections, which shift the scalar glueball mass by about few percent [20], will be neglected in this work. Other theoretical approaches to the exploration of glueball physics have been introduced in the review [21]. Note that there exist several estimates for a crucial higher-power contribution to the OPE, i.e., the triple-gluon condensate, such as those from the single-instanton [1,22,23] and lattice [24] evaluations, and that determined from heavy quark systems [25], which differ dramatically. Once the spectral density is obtained from the dispersion relation, one can compute the correlation function associated with the scalar glueball at zero momentum. It will be elaborated that the spectral density corresponding to the lattice estimate respects better the low-energy theorem for the correlation function [14,26,27]. That is, our method helps discriminate the different estimates for the triple-gluon condensate.
Appropriate moments of a Borel transformed correlator have to be selected to form ratios for the extraction of a glueball mass in sum rules, because the stability window in the Borel mass may not exist for ratios of other moments [10,13]. On the other hand, the lower (higher) moments are more sensitive to light (heavy) resonances [13]. The above subtlety is not an issue to our formalism, in which full information of the OPE input is utilized to solve a dispersion relation. It will be seen that a double-peak structure with both light and heavy resonances is revealed by the solution to a spectral density, in support of the observation in [28]. We emphasize that the positivity constraint on the spectral density plays a crucial role for establishing the bound states. The shorter peak located at 0.60 GeV in the spectral density for the scalar glueball can be identified as the light scalar meson f 0 (500) perhaps with little admixture from f 0 (980). Another located at 1.50 GeV with a broader width about 200 MeV arises from the combined contribution of the f 0 (1370), f 0 (1500) and f 0 (1700) mesons, since f 0 (1500) with a narrow width cannot accommodate the broad peak alone. Our solution indicates that f 0 (500) and f 0 (980) contain some amount of gluonia, and the f 0 (1370), f 0 (1500) and f 0 (1700) mesons are the glue-rich states, in agreement with the multiple-resonance mixing picture for scalar mesons widely discussed in the literature [29][30][31][32][33][34][35]. The prediction for the scalar glueball mass around 1.50 GeV is close to those from sum rules [13,18,36] (but a bit lower than in [37]) and quenched lattice QCD [38][39][40][41], and smaller than from holographic QCD [42]. The quenched lattice calculation in [43] favors the interpretation of f 0 (1710) as composed mainly of the lightest scalar glueball. We mention that a double-peak parametrization for the spectral density has been found to yield a fit to the sum rule for the scalar glueball with the instanton effect better than a single-peak one [11].
The spectral density for the pseudoscalar glueball is also featured with a double-peak structure: the shorter peak located at 0.71 GeV comes from the combined contribution of the η and η mesons naturally, which have been known to comprise some gluonium components [44]. Another at 1.75 GeV with a broad width about 200 MeV strongly suggests that the η(1760) meson is a promising candidate for the pseudosclar glueball. This mass is lower than most results in the literature, which are above 2 GeV, such as those from sum rules [13], quenched lattice QCD [38][39][40][41], the Bethe-Salpeter approach [45,46] and holographic QCD [42]. Nevertheless, when the resonance contribution was parametrized by a Breit-Wigner form with a finite width, the pseudoscalar glueball mass drops to 1.407 ± 0.162 GeV in sum rules with the instanton effect [19]. We remind that measurements of J/ψ radiative decays do not affirm any glue-rich pseudoscalar resonances with masses above 2 GeV (the quantum numbers of X(2370) are not certain yet) [12]. The η(1760) meson was proposed to be the pseudoscalar glueball two decades ago [47,48], and examined experimentally via the decay J/ψ → γ(η(1760) →)ωω in [49]. It is abundantly produced in J/ψ radiative decays, but not seen in the J/ψ → γγV channels [50], V = ρ, φ, implying that the partial width of η(1760) → γV is tiny, namely, η(1760) is likely to be a glueball. Our solution to the spectral density for the pseudoscalar glueball advocates that η(1760) should be included into the mixing framework for pseudoscalar mesons [51][52][53][54][55] for a complete scenario.
The topological susceptability χ t , an important quantity characterizing the nonperturbative QCD vacuum, is related to the low-energy limit of the correlation function for the pseudoscalar glueball, which has been analyzed in sum rules with the pure Yang-Mills theory [56]. Similarly, once the spectral density is ready, one can compute the corresponding correlation function at zero momentum. We then predict χ 1/4 t = 75-78 MeV, which is quite precise and matches those from lattice QCD [57][58][59][60][61][62] and the chiral perturbatin theory [63][64][65]. Viewing that our formalism produces the observed ρ meson mass, the widely accepted scalar glueball mass around 1.5 GeV, and the topological susceptability in consistency with the expectations from other approaches, we are convinced that the η(1760) meson deserves thorough investigations on whether it is the pseudoscalar glueball in the long-lasting quest.
At last, we find no resonance solution to the spectral density for the tensor glueball. Distinct from the scalar and pseudoscalar glueball cases, only a single dimension-eight condensate is available for the OPE input. It has been verified [2] that at least two condensates of different dimensions are necessary for establishing the ρ meson state. We thus speculate that the absence of a resonance solution for the tensor glueball may be due to the insufficient nonperturbative input in the present setup. Nevertheless, a tensor glueball mass about 2.0 GeV was extracted from specific moments of the sum rule in [13]. As stated before, full information of the OPE input is utilized in our formalism, so the condition on the existence of a resonance may be more stringent. It is urged to calculate more higher-power corrections to the OPE of the correlation function, so that the tensor glueball mass can be inferred from the dispersion relation and compared with the lattice predictions [39,66,67].
The rest of the paper is organized as follows. In Sec. II we illustrate the inverse matrix method to solve a dispersion relation as an inverse problem. The mock data from several test functions are constructed, and treated as inputs. It is shown that the solutions reproduce the test functions accurately in the inverse matrix method. Our approach is then applied to the determination of the ρ meson mass, via which the value of the gluon condensate is fixed for the analyses of the glueball masses. We explain how boundary conditions of a spectral density select the suitable set of generalized Laguerre polynomials for its expansion, and why it is difficult to probe excited states. In Sec. III we solve for the scalar, pseudoscalar and tensor glueball masses, suggest the physical states they correspond to, and discuss their impact on meson mixing scenarios. The correlation functions at zero momentum for the scalar and pseudoscalar glueballs are also deduced from the spectral densities. The former serves to discriminate the different estimates for the triple-gluon condensate, and the latter is used to predict the topological susceptability. Section IV contains the conclusion and outlook.

A. Inverse Matrix Method
We first demonstrate our inverse matrix method to solve a dispersion relation, which belongs to the first kind of Fredholm integral equations, by means of several simple examples. The goal is to find the unknown function ρ(y) from the integral equation given the input function ω(x) of the variable x. Suppose that ρ(y) decreases quickly enough with the variable y, so the major contribution to the integral on the left-hand side originates from a finite range of y. It is then justified to expand the integral into a series in 1/x up to some power N for a sufficiently large |x| by inserting into Eq. (1). Also suppose that ω(x) can be expanded into a power series in 1/x for a large |x|, Note that |x| being large enough is only a formal statement, and does not have a substantial influence on our calculation.
There are four types of classical orthogonal polynomials constructed from solutions to second-order differential equations: the Jacobi-like polynomials (including the Gegenbauer polynomials, the Legendre polynomials,...) with the support [−1, 1]; the Laguerre polynomials with the support [0, ∞), the Bessel polynomials with the support [0, ∞), and the Hermite polynomials with the support (−∞, ∞) [68]. Since ρ(y) in Eq. (1) (and also spectral densities to be studied later) is defined in the interval [0, ∞), the Laguerre and Bessel polynomials may be considered for its expansion. However, the latter are orthogonal with respect to the weight exp(−2/y), which suppresses the small-y region we are interested in. Hence, the Laguerre polynomials, orthogonal with respect to the weight exp(−y), are the only appropriate choice for the basis functions. We thus decompose the unknown into in terms of a set of generalized Laguerre functions L (α) n up to degree N − 1, which satisfies the orthogonality The maximal number N will be fixed later, and the choice of the parameter α depends on the behavior of ρ(y) in the limit y → 0. Substituting Eqs.
with m and n running from 1 to N , and the vectors a = (a 1 , a 2 , · · · , a N ) and b = (b 1 , b 2 , · · · , b N ). If the inverse of M exists, one can get a solution of a via a = M −1 b with the known input b trivially. The existence of M −1 thus implies the uniqueness of the solution to ρ(y). In principle, the true solution can be approached to by increasing the number of polynomials N in Eq. (4). The difference between the true solution and the approximate one produces a power correction 1/x N +1 to the left-hand side of Eq. (1), because of the orthogonality in Eq. (5), which is beyond the considered precision. The orthogonality also leads to M mn = 0 for m < n. Namely, M is a triangular matrix, such that the coefficients a n built up previously are not altered, when an additional higher-degree polynomial is added to the expansion in Eq. (4). Nevertheless, both m and n have to stop at a finite N in a practical application, since the determinant of M diminishes with its dimension eventually. An approximate solution of a would then deviate from the true solution violently, when a tiny fluctuation of the input vector b is present and amplified by the huge elements of M −1 . This is a generic feature of an ill-posed inverse problem. Hence, the optimal choice of N is set to its maximal value, above which a solution goes out of control.
We test the above inverse matrix method on the following simple examples. We generate the mock data using a single-peak function for the input ω(x) with the coefficients The factor e −y 2 in the test function guarantees that the dominant contribution to the integral comes from the region with finite y, and justifies the power expansion in 1/x with a sufficiently large |x|. The matrix M is computed according to Eq. (6) with α = 2, motivated by the limit ρ(y) → y 2 as y → 0. The inverse M −1 gives the solutions for the coefficients a n via a = M −1 b, and ρ(y) in Eq. (4), whose behaviors with the expansions up to N = 22, 23 and 24 generalized Laguerre polynomials L  23 are almost identical, implying the stability of the solutions for a sufficiently large N . In fact, the solutions have changed little as N > 20. It is also seen that the curve becomes oscillatory, and differs drastically from the single-peak function when N reaches 24. The big ratio of the last two coefficients a 24 /a 23 ≈ 58, compared with a 22 /a 21 ≈ 1 for N = 22 and a 23 /a 22 ≈ 2 for N = 23, hints that the inverse M −1 is out of control at N = 24.
We fix N = 22 in the demonstration below, though the choices N = 21 and N = 23 also serve the purpose. Figure 2 collects the solutions to ρ(y) from the parameters α = 0, 2 and 4 for the generalized Laguerre polynomials L (α) n (y), compared with the true solution in Eq. (7). The curve corresponding to α = 2, namely, the expansion of ρ(y) in Eq. (4) with the correct behavior in the y → 0 limit, matches the true solution most perfectly. The curve labelled by α = 0, despite of describing the true solution equally well at finite y, shows deviation near the origin y = 0. As α = 4, the approximate solution differs completely from the true solution. The above test manifests the importance of the information on the boundary conditions of the unknown and the choice of the appropriate set of generalized Laguerre polynomials for solving Eq. (1). It will be made explicit in the next sections that the α = 1 (α = 2) set of generalized Laguerre polynomials is selected for the analysis of the ρ meson (glueball) mass.
We highlight that the exponential factor e −y in the polynomial expansion in Eq. (4) characterizes the resolution power ∆y ∼ 1 of our method to probe the structure of the unknown ρ(y). To elaborate this point, consider the two double-peak functions in which the two peaks are separated by ∆y ∼ 0.5 and ∆y ∼ 1.5, respectively. The large coefficients 20 in the exponents are designed to make sharp peaks for transparent illustration. The above two functions are substituted for y 2 e −y 2 in Eq. (8) to generate the mock data for the inputs. The same procedure, with the expansion up to 22 generalized Laguerre polynomials L (0) n (y), gives the solutions to ρ 1 (y) and ρ 2 (y) in Fig. 3. It is observed that the approximate solution contains a single peak in the left plot, since the two peaks of ρ 1 (y) are too close to resolve. On the contrary, the two-peak structure of ρ 2 (y), with a separation within the designated resolution, can be reproduced reasonably. The above test reflects the limitation of our method to probe a physical system with certain characteristic scales.

B. ρ Meson Mass
After exploring the aspects of our approach, we apply it to the determination of the ρ meson mass from the corresponding dispersion relation. We first recapitulate and expand the idea of handling QCD sum rules as an inverse problem [2], starting with the two-point correlator for the quark current J µ = (ūγ µ u −dγ µ d)/ √ 2. The vacuum polarization function Π(q 2 ) obeys the identity where the contour consists of two pieces of horizontal lines above and below the branch cut along the positive real axis on the complex s plane, and a circle of large radius R [2]. The OPE of the function Π(q 2 ) in the deep Euclidean region of q 2 is reliabe, and we have Π OPE (q 2 ) [1] for the left-hand side of Eq. (11), up to the dimension-six condensate, ie., to the power correction of 1/(q 2 ) 3 . The second expression of the perturbative term Π pert (q 2 ) defines the constant c. In Eq. (12) α s G 2 ≡ α s G a µν G aµν is the gluon condensate, m q is a quark mass, and the parameter κ = 2-4 [69][70][71] quantifies the violation in the factorization of the four-quark condensate (qq) 2 into the square of the quark condensate qq . A regularization-scheme dependent constant in Π pert (q 2 ) [72] has been dropped, which is irrelevant to the search for a resonance solution. The fact that this constant can be eliminated by the Borel operator in standard sum rules confirms the above statement.
The contour integral on the right-hand side of Eq. (11) can be written as in which the lower bound of the first integral on the right-hand side, being of order of the pion mass squared, has been set to zero for simplicity, and the imaginary part ImΠ(s), involving nonperturbative dynamics from the low s region, will be solved for later. The numerator in the second integral, with C representing the large circle of radius R, has been replaced by Π pert (s), because the perturbative evaluation of Π(s) is reliable for s far away from physical poles, in accordance with the employment of the OPE in Eq. (12). We also express the perturbative piece Π pert (q 2 ) by means of an integration along the same contour, so Eq. (12) becomes Equating Eqs. (14) and (15) according to Eq. (11), we get the sum rule where the contributions of Π pert (s) along the big circle C, together with the dependence on the renormalization scale µ, have cancelled from both sides. We introduce a subtracted spectral density, which is related to the original one ρ(s) ≡ ImΠ(s)/π via The scale Λ characterizes the transition of ImΠ(s) to the perturbative expression ImΠ pert (s). The smooth function 1 − exp(−s/Λ) vanishes like s as s → 0, and approaches to the unity at large s Λ, such that ∆ρ(s, Λ) respects the behavior of ρ(s) ∼ s in the limit s → 0 [73], and diminishes quickly as s > Λ. Note that ∆ρ(s, Λ) bears the nontrivial resonance structure the same as ρ(s) for s < Λ, which is not affected by the perturbative subtraction term. We have confirmed that other smooth functions with the similar limiting behaviors lead to basically identical solutions for ρ(s). If one adopts the step function, instead of the smooth function in Eq. (17), to define ∆ρ(s, Λ), the resultant ρ(s), as a sum of the smooth ∆ρ(s, Λ) from solving the dispersion relation and the discontinuous perturbative contribution caused by the step function, will exhibit a sudden jump. The radius R in Eq. (16) can be pushed toward the infinity, when the sum rule is formulated in terms of the subtracted spectral density: where the constant c has been defined in Eq. (13), and the aforementioned regularization-scheme dependent constant is absent from ImΠ pert (s). It is stressed that the quark-hadron duality for the unknown spectral density is not assumed at any finite scale s in the above derivation.
Since the subtracted spectral density ∆ρ(s, Λ) is a dimensionless quantity, it can be expressed as a function in the form ∆ρ(s/Λ). Certainly, the subtracted spectral density may depend on other constant scales, such as the ρ meson mass m ρ , which, however, appears as a constant ratio m ρ /Λ for a given Λ, and needs not be shown as an explicit argument. Equation (18) then reduces, with the variable changes x = q 2 /Λ and y = s/Λ, to where Λ in the subtracted spectral density has moved into the condensate terms to make them dimensionless. On one hand, the scale Λ prefers to be small to enhance the resolution of our method. On the other hand, it cannot be too small to spoil the OPE. When Λ increases from a low scale, the physical ρ meson mass, if generated, corresponds to a peak location of ρ(s), which should be insensitive to the change of Λ. As Λ further increases, it disappears with the condensate contributions from Eq. (19). Then a solution for ∆ρ(y), if existent, will imply that ∆ρ(s/Λ) is a solution for an arbitrary Λ. A peak location of ρ(s), endowed by ∆ρ(s/Λ), thus shifts with Λ, such that none of its structure can be interpreted as a physical state. The numerical analyses to be performed below verify this tendency of the ρ meson mass m ρ , as well as of the glueball masses, with respect to the variation of Λ: the obtained m ρ remains stable first, and then grows with Λ monotonically. In the sense of searching for a stability window in which a resonance mass stays constant, Λ plays a role similar to the Borel mass in conventional sum rules.
Viewing the boundary conditions of ∆ρ(y) ∼ y at y → 0 and ∆ρ(y) → 0 at y → ∞, we expand ∆ρ(y) in terms of the generalized Laguerre polynomials L (1) n (y), namely, employ Eq. (4) with the parameter α = 1. The dependence on the constant ratios mentioned before then goes into the coefficients a n : a solution of a n depends on Λ, as indicated by the right-hand side of Eq. (19). Equation (2), which holds for |x| > 1, ie., for |q 2 | > Λ, is inserted into the left-hand side of Eq. (19) to construct the matrix elements M mn in Eq. (6), and inserted into the integral on the right-hand side of Eq. (19) to compute the coefficients b n for the input. The coefficients b 2 and b 3 of the 1/x 2 and 1/x 3 terms, respectively, receive additional contributions from the condensates. The following OPE parameters and the strong coupling α s , evaluated at the scale of 1 GeV and within their accepted ranges [1, 13, 16, 69-71, 74, 75], are adopted: which will be shown to produce the observed ρ meson mass m ρ = 0.78 GeV [12]. Though higher-order corrections to the condensate terms are available [74,76], it may not mean much to include them, because their effects can be mimicked by tuning the unknown factorization violation parameter κ. We have checked that the renormalizationgroup evolutions of α s and the condensates around the scales 1-2 GeV affect our results for m ρ by only few percent, so they will not be taken into account in the numerical study. We derive the inverse matrix M −1 , the coefficients a n from the OPE coefficients b n , the solution to ∆ρ(s, Λ), and then the spectral density ρ(s) from Eq. (17). The outcomes for the characteristic scale Λ = 2.5 GeV 2 with the expansions up to N = 22, 23, and 24 generalized Laguerre polynomials L (1) n (y) are displayed in Fig. 4. It is found that the curves labelled by N = 22 and 23 are very similar, assuring the stability of the solutions for a sufficiently large N , and consistent with those obtained in the the maximum entropy method [9]. The curve becomes oscillatory, and differs significantly from the other two as N reaches 24. The big ratio of the last two coefficients a 24 /a 23 ≈ 7, compared with a 22 /a 21 ≈ a 23 /a 22 < 2 for N = 22 and 23, indicates that the matrix elements of M −1 start to increase rapidly as N = 24. It is encouraging that the positivity of the spectral density is satisfied automatically. We read the ρ meson mass m ρ = 0.78 GeV (m 2 ρ = 0.61 GeV 2 ) off the location of the sharp peak in the plot for N = 23, which agrees with the measured value in [12]. The bump located at s > 2 GeV 2 , being shorter and broader, may be attributed to the combination of nonresonant contributions and resonant ones from excited ρ states. It is obvious that the continuum contribution to the spectral density ρ(s), distinct from the perturbative value c ≈ 0.029, violates the local quark-hadron duality, though ρ(s) approaches to c asymptotically. We then seek the solutions of Eq. (19) for various characteristic scales Λ = 1-8 GeV 2 by repeating the above procedure, and present the spectral densities ρ(s) for Λ = 1, 4 and 8 GeV 2 in Fig. 5 for illustration. It is noticed that the expansion in terms of the generalized Laguerre polynomials diverges more quickly at lower Λ, such that solutions for Λ < 1 GeV 2 with the polynomial numbers roughly smaller than 10 may not be reliable. A clear signal for a divergent expansion is that the positivity of the spectral density is lost owing to strong oscillation. We have to terminate the expansion at N = 12 (N = 22) for Λ = 1 GeV 2 (Λ = 2 GeV 2 ), and at N = 23 for higher Λ. Figure 5 suggests that the peaks are located at the same s ≈ 0.61 GeV 2 for Λ between 1 and 4 GeV 2 , which specifies the stability window in Λ for the physical solutions. The peak in the plot for Λ = 8 GeV 2 moves toward a bigger s = 1.24 GeV 2 , exhibiting the growth of the mass at large Λ with the disappearance of the condensate effects. We point out that the peak of ρ(s) always exists even without the condensates [2]: the vanishing of the spectral density at s = 0 gives a contribution to the dispersive integral smaller than the constant perturbative piece does, which must be made up by a peak at finite s in order to match the OPE input. Therefore, one has to be cautious about the identification of a peak in ρ(s) as a physical resonance, for which the stability of the peak location inspected above is crucial. Accordingly, the broad bump in Fig. 5 cannot be interpreted as a specific bound state, since its location shifts with Λ as disclosed in the three plots. We also observe that the peak becomes less evident at large Λ, for the bad resolution of our method stops probing the structure of ρ(s) at low s. Next we scan the ρ meson mass m ρ in the range 1 GeV 2 < Λ < 8 GeV 2 , and depict its dependence on Λ in Fig. 6. The expected features are salient: the curve goes up and down around m ρ = 0.78 GeV in the interval 1 GeV 2 < Λ < 3 GeV 2 , and then ascends monotonically with Λ as Λ > 3 GeV 2 . We assess the theoretical error from our method using the extreme values at Λ = 2.2 GeV 2 and 3.0 GeV 2 , and get m ρ = (0.78 ± 0.03) GeV. To check the sensitivity of our results to the uncertainties of the OPE, we vary the gluon condensate α s G 2 and the factorization violation parameter κ by ±20% separately. The variation of the former can simulate that of the quark condensate qq , which is also of dimension four. The variation of the latter is equivalent to that of the dimension-six condensate qq 2 . It turns out that the obtained m ρ changes by only about ∓5% and ±7%, respectively. Namely, results in our formalism are insensitive to the OPE uncertainties. Adding the above sources of errors in quadrature, we conclude m ρ = (0.78 ± 0.07) GeV.
The above examination also reveals that α s G 2 and κ must be anti-correlated in order to fix m ρ . Hence, the input value of κ ( α s G 2 ) in Eq. (20), close to the lower (upper) edge of its favored range, is the optimal choice. It was stated [77] that the area under the resonance peak of the spectral density ρ(s) represents the square of the ρ meson decay constant f ρ . In our formulation, the resonance peak is appropriately described by the subtracted spectral density ∆ρ(s, Λ), where the continuum contribution has been largely removed. We thus have with the scale Λ = 1.3 GeV 2 corresponding to the ρ meson mass m ρ = 0.78 GeV in Fig. 6. The above integral leads to f ρ ≈ 0.20 GeV, consistent with the value in [12], which further supports our approach to the extraction of nonperturbative observables from dispersion relations. Note that f ρ derived from Eq. (22) depends on Λ: varying Λ from 1.0 to 3.0 GeV 2 , we find that f ρ increases from 0.17 to 0.29 GeV. This Λ dependence of f ρ is expected, because Eq. (22) holds better for a narrower resonance, but the ρ meson width is not small: Eq. (22) can be understood via a pole parametrization, which originates from the narrow-width approximation. At last, a remark is in order. The well-known Weinberg sum rules [78] for the difference ρ V (s) − ρ A (s), where ρ V (s) corresponds to the vector spectral density investigated in this subsection and ρ A (s) denotes the axial-vector spectral density, have been analyzed in the literature. These sum rules, according to their derivation, are expected to be respected in our formalism, where the spectral densities are solved with the input from the OPE of the relevant correlation functions. Certainly, one has to obtain ρ A (s) before verifying the above postulation, which can be a subject of future studies. In such an approach, the Weinberg sum rules serve as a check of the theoretical consistency of our solutions for the spectral densities. It differs from the approaches in [79,80], where the data of e + e − annihilation and τ decays were employed to construct ρ V (s) and ρ A (s) for examining whether the sum rules are satisfied.

A. Formalism
We apply the formalism developed in the previous section to the determination of the scalar and pseudoscalar glueball masses. Consider the correlation function for the glueball channel where the local composite operators O G with G = S and P denote the gluonic interpolating fields for the scalar (0 ++ ) and pseudoscalar (0 −+ ) glueballs, respectively. Their explicit definitions with the lowest mass dimension are given by G µν ≡ i µνρσ G ρσ /2 being the dual of the gluon field strength. The low-energy theorem demands the zero-momentum limits of the glueball correlators [14,26]: in which β 0 = 11N c /3 − 2N f /3, N c and N f being the numbers of colors and flavors, respectively, is the lowest-order coefficient of the QCD β-function, and the topological susceptibility is defined with the topological charge density Q(x) = G aµν (x)G a µν (x)/(32π 2 ). The above low-energy limits will be used to test the consistency of our calculations below.
Similarly, we have the OPE of the correlation function Π G (q 2 ) in the deep Euclidean region of q 2 [1], up to the dimension-eight condensate, ie., up to the power correction of 1/(q 2 ) 2 , where the various gluon condensates are defined as The four-gluon condensates are approximated, under the vacuum factorization assumption [81,82] ( by As stated in the introduction, we consider only the condensate contributions to the OPE, which are sufficient for establishing the scalar and pseudoscalar glueballs, without the instanton effect. The coefficients in Eq.
Those for the pseudoscalar glueball were found to be [81,87,88] Notice C (S) 1 = C (P ) 1 = 0, which will not appear in our formulas afterwards.
We extend the derivation of Eqs. (11)- (16) to the construction of the dispersion relations for the glueball masses, arriving at The imaginary part ImΠ pert G (s) collects the contributions in Eq. (28) without poles at q 2 → 0, ie., those which can be produced by the contour integration of the perturbative piece, like the first term in Eq. (15): It is seen that the term B (G) 0 α s G 2 in Eq. (28) is absent in the above expression, since it has no discontinuity along the branch cut. As stated before, a constant piece in the OPE is irrelevant to the search for a resonance solution. The thresholds in the dispersive integrals on both sides of Eq. (34) have been set to zero for simplicity. We remind that ImΠ pert G (s) contains the nonperturbatve gluon condensate α s G 2 actually. We introduce the subtracted spectral density in terms of the original one ρ G (s) ≡ ImΠ G (s)/π. The smooth function 1 − exp(−s/Λ) vanishes like s and 1 − exp(−s 2 /Λ 2 ) vanishes like s 2 as s → 0, so the subtraction terms do not modify the behavior of ρ G (s) ∼ s 2 in the low-energy limit s → 0 [14,26,89,90]. Note that s 2 ln s and s 2 ln 2 s are regarded as being much larger than s 2 as s → 0, so additional suppression from the first smooth function is necessary. Both the smooth functions approach to the unity at large s Λ, where the subtracted spectral density ∆ρ(s, Λ) diminishes quickly, and the radius R can be pushed toward the infinity. The scale in the exponent of the second smooth function can differ from Λ in principle. However, we have verified that our numerical results are insensitive to its variation, so it is chosen as Λ without loss of generality. It has been also confirmed that the solutions to ρ G are basically identical, when different smooth functions are employed in Eq. (36). The subtracted spectral density ∆ρ G (s, Λ) exhibits the resonant structure almost the same as ρ(s) does, which is barely affected by the subtraction term. It will be shown that the introduction of a subtracted spectral density facilitates the evaluation of the correlation function at zero momentum.

Equation (34) is then converted into
The renormalization scale µ is usually set to the Borel mass, when the Borel transformation is applied to sum rules [10,13,16]. We set µ 2 to the characteristic scale Λ, and apply the variable changes x = q 2 /Λ and y = s/Λ to Eq. (37), obtaining where the dimensionless function ∆ρ G (s, Λ)/Λ 2 has been replaced by ∆ρ G (y), according to the reasoning in the previous section. The scale Λ characterizes the region with y < 1, from which the dominant nonperturbative contribution to the dispersive integral arises. It has been argued that a range of Λ exists, in which a glueball mass m G , corresponding to a peak location of ∆ρ G (s, Λ), is stable against the variation of Λ. As Λ becomes large enough, it diminishes the nonpertrbative condensate effects, and the scaling of a solution with Λ appears. When the scaling occurs, no structure of a solution can be interpreted as a physical state. Since the subtracted spectral density ∆ρ G (y) follows the behaviors ∆ρ G (y) ∼ y 2 as y → 0 and ∆ρ G (y) → 0 as y → ∞, we expand it in terms of the generalized Laguerre polynomials L  Equation (2) is inserted into the left-hand side of Eq. (38) to compute the matrix elements M mn in Eq. (6), and inserted into the right-hand side of Eq. (38) to gain the coefficients b n for the input. The coefficients b 1 and b 2 of the 1/x and 1/x 2 terms, respectively, receive additional contributions from the condensates. We take the inputs of the gluon condensate α s G 2 and of the strong coupling α s the same as in Eq. (20) for consistency. The triple-gluon condensate is given by from the single-instanton estimate [1,22,23] and the lattice estimate [24], respectively. The value in [25] has a positive sign with magnitude about five times higher than in Eq. (39). They are different apparently as having been noticed in [91], and affect numerical outcomes. We will discriminate the above estimates by the low-energy limit of the correlation function for the scalar glueball in Eq. (25). Besides, we do not consider the minor renormalization-group effect either. As emphasized before, the low-energy limit of the correlation function can be derived from a solution to the spectral density. Start with the dispersion relation for the correlation function Π G (q 2 ) similar to Eq. (14), and write in which Eq. (36) has been substituted for ImΠ G (s)/π, and the pieces in Eq. (35) have been grouped into the second term on the right-hand side to allow a closed contour. Because the dispersive integral is finite, the radius R can be pushed to infinity. We then derive the low-energy limit where the last term proportional to α s G 2 is from the contour integral for q 2 = 0 and µ 2 = Λ. The constant B (G) 0 , despite of being absent in Eq. (38), appears in the above expression. We point out that ∆ρ G (y) corresponds exactly to the ultraviolet regularized spectral density with the high-frequency contribution being removed, which has been used to define a correlation function at zero momentum in [16]. Equation (42) thus specifies explicitly how to perform the ultraviolet regularization for a dispersive integral.

B. Scalar Glueball Masses
We extract the scalar glueball mass by solving Eq. (38) as an inverse problem first with the triple-gluon condensate in Eq. (39), getting the inverse matrix M −1 , the coefficients a n from the OPE coefficients b n , the solution to ∆ρ S (s, Λ), and the spectral density ρ S (s) for various characteristic scales Λ. The results of ρ S (s) for Λ = 1.5 GeV 2 with the n (y) are displayed in Fig. 7. The coefficients a n have not yet grown quickly, and the three curves are similar, implying the stability of the solutions. We find the ratios of the last two coefficients, a 15 /a 14 ≈ a 16 /a 15 ≈ a 17 /a 16 ≈ 1 in the three cases. However, the spectral density ρ S (s), supposed to be positive, becomes negative around s ≈ 1.2 GeV 2 for N = 17. The positivity constraint forces us to terminate the polynomial expansion at N = 16. That is, the positivity of the spectral density plays a more important role in the determination of glueball masses than in the ρ meson case. The same procedure for the triple-gluon condensate in Eq. (40) yields the solutions for the characteristic scale Λ = 1.5 GeV 2 with the expansions up to N = 18, 19 and 20 generalized Laguerre polynomials L (2) n (y) in Fig. 8. The curves in the three plots are also similar, and we select the one with N = 19 as the solution to avoid violating the positivity constraint.
A closer look reveals that the two solutions with N = 16 in Fig. 7 and with N = 19 in Fig. 8 are different: for instance, the first peak of the former is taller and located at larger s ≈ 0.5 GeV 2 , while the first peak of the latter is shorter and located at lower s ≈ 0.4 GeV 2 . It means that one has to decide which triple-gluon condensate is adopted in order to predict glueball masses unambiguously. The discrimination can be achieved by confronting the solutions with Eq. (25) from the low-energy theorem, namely, Π S (0) = 0.78 GeV 4 for the gluon condensate α s G 2 = 0.08 GeV 4 in Eq. (20). Inserting the subtracted spectral densities ∆ρ S (s) corresponding to Eqs. (39) and (40) into Eq. (42), we have Π S (0) = 1.08 and 0.60 GeV 4 , respectively. That is, the lattice estimate for the triple-gluon condensate in Eq. (40) leads to the zero-momentum correlation function more consistent with the low-energy theorem. If the large triple-gluon condensate in [25] is adopted, we will find Π S (0) ≈ 1.7 GeV 4 for Λ = 1.5 GeV 2 , and that the first peak of the spectral density shifts to a higher s ≈ 0.6 GeV 2 . The central value of Π S (0) resulting from the same triple-gluon condensate was also found to be larger than that set by the low-energy theorem recently [36]. It has been ensured by varying the characteristic scale Λ that the values Π S (0) from the single-instanton estimate in Eq. (39) and from [25] are far above 0.78 GeV 4 (Π S (0) is not sensitive to the variation of Λ). Therefore, we will pick up the input in Eq. (40) for the rest of numerical investigations.
It is interesting to notice in Fig. 8 that two peaks, one at s ≈ 0.4 GeV 2 and another at s ≈ 2.3 GeV 2 , appear in all the three plots, though the latter seems not obvious due to the huge perturbative background from ImΠ perp S (s)/π. To highlight these two peaks, we turn to the subtracted spectral density ∆ρ S (s, Λ) below. The expansion in terms of the generalized Laguerre polynomials diverges more quickly at lower Λ, such that solutions for Λ < 1 GeV 2 with the polynomial numbers roughly smaller than 10 may not be reliable. The maximal number N for the polynomial expansion then increases with Λ under the positivity constraint till Λ = 1.6 GeV 2 , above which the spectral density is always positive, but the matrix elements of M −1 and the coefficients a n grow quickly. For example, we find a 23 /a 22 ≈ 2.3 for Λ = 2.5 GeV 2 and N = 23, so we stick to N = 22 for the polynomial expansion in the range Λ = 1.7-2.5 GeV 2 . The above prescribes how the maximal number of polynomials is fixed for an expansion. The behaviors of ∆ρ S (s, Λ) for Λ = 1.5, 1.8 and 2.5 GeV 2 are shown in Fig. 9, where the double-peak structure in the first two plots is significant. We stress that these two peaks are not a numerical artifact, since we have demonstrated that a double-peak structure can be disclosed by our method, as the resolution power, ie., the scale Λ is appropriately chosen. In the present case Λ ∼ O(1) GeV 2 is able to resolve the two peaks located at m 2 S1 ≈ 0.4 GeV 2 and m 2 S2 ≈ 2.3 GeV 2 . The similar shapes in the first two plots of Fig. 9 for Λ = 1.5 and 1.8 GeV 2 indicate that the solutions for the scalar masses are stable in an interval of Λ. Note that the smooth functions in the subtraction term in Eq. (36) are close to the unity as s > 4 GeV 2 for the considered Λ = 1.5 or 1.8 GeV 2 . Therefore, the small deviation from zero above s ≈ 4 GeV 2 signals the mild violation of the local quark-hadron duality of the spectral density ρ S (s) in the high s region. It is apparent that the third plot in Fig. 9 labelled by Λ = 2.5 GeV 2 differs much from the first two: the peaks move toward larger s, manifesting the scaling behavior of ∆ρ S (s, Λ) at large Λ ascribed to the disappearance of the nonperturbative condensate effects. The double-peak structure also blurs, since the resolution of our method becomes worse at this large Λ.
We remark that the matrix element M mn in Eq. (6) corresponds to the m-th moment of a sum rule for a glueball mass. As postulated in [13], the lower moments of a sum rule are more sensitive to low-lying resonances, while higher moments are, on the contrary, more sensitive to heavy resonances: the light masses about 700-900 MeV were extracted from the lower moments of the sum rule in [92], and the heavy ones about 1.5-1.7 GeV were extracted from the higher moments in [10,13,84]. It is then easy to understand why our formalism, which takes into account more moments than in conventional sum rules [13,93,94], produces the two peaks simultaneously with the similar squared masses m 2 S1 ≈ 0.4 GeV 2 and m 2 S2 ≈ 2.3 GeV 2 . The double-peak structure in Fig. 9 is consistent with the observation in the sum-rule analysis [11], where a double-resonance parametrization for the spectral density was shown to give a fit to the OPE side better than a single-resonance one. The mass range 0.8-1.6 GeV for the peak locations in [11] is similar to ours, but the mass gap between the two resonances, about few hundreds of MeV in [11], is smaller than in our solution, which is about 1 GeV. Another distinction is that the lighter resonance has a broader width in [11], while the heavier one does in our solution, which may be due to the different nonperturbative effects included in the theoretical frameworks: they come from the instanton contribution to the correlation function in [11], but from the gluon condensates in ours. We mention that the scalar glueball mass around 1.5 GeV, very close to our m S2 , was derived from sum rules [15,18] with only the instanton effects. It supports our claim that the gluon condensates are sufficient for establishing the glueballs.
We scan the scalar masses m S1 and m S2 in the range 1 GeV 2 < Λ < 2.5 GeV 2 , and depict their dependencies on Λ in Fig. 10, where the two curves describe the two peak locations of the subtracted spectral density ∆ρ S (s, Λ). It is seen that the lower and upper curves ascend first from Λ = 1 GeV 2 , reaching m S1 = 0.60 GeV and m S2 = 1.50 GeV, respectively, become stable in the window 1.2 GeV 2 < Λ < 1.7 GeV 2 , and then go up again monotonically. These features, having appeared in Fig. 6 for the ρ meson mass, are completely anticipated. We then survey the dependence of the correlation function Π S (0) at zero momentum on Λ according to Eq. (42), given the solutions for ∆ρ S (s, Λ), and obtain Π S (0) = 0.65 GeV 4 at Λ = 1.7 GeV 2 . Considering the typical 20% uncertainty around Π S (0) = 0.78 GeV 4 from Eq. (26), we are sure that the low-energy limit is satisfactorily respected by our solutions in the stability window. It is legitimate to choose m S1 = 0.60 GeV and m S2 = 1.50 GeV as the central values, and to estimate the theoretical errors in our method using the minimal (maximal) values at Λ = 1.4 (1.5) GeV 2 . We get m S1 = (0.60±0.01) GeV and m S2 = (1.50 ± 0.01) GeV, whose tiny errors reflect the remarkable stability of our solutions. Because the subtracted spectral densities for other values of Λ in the interval 1.2 GeV 2 < Λ < 1.7 GeV 2 are very similar to that for Λ = 1.5 GeV 2 in Fig. 9, we do not present them here.
We investigate the theoretical uncertainties arising from the variation of the involved parameters, which all turn out to be under control. Decreasing the strong coupling α s from 0.5 to 0.4, which corresponds to the scale variation within the stability window roughly, we find 8% enhancement on the determined scalar masses. This check justifies the neglect of the renormalzation-group effect in our calculation. The typical ±20% change of the gluon condensate α s G 2 causes about ±5% impact. We also examine the sensitivity of the scalar masses to choices of the renormalization △ρ P /Λ 2 s (GeV ) scale µ, and observe 2% increase from µ 2 = 2Λ and µ 2 = Λ/2. Adding the above sources of errors in quadrature, we conclude m S1 = (0.60 ± 0.06) GeV, m S2 = (1.50 ± 0.15) GeV.
We interpret the solutions to the spectral density for the scalar glueball, bearing in mind that any physical state with a gluonic content can contribute to the considered spectral density. The major peak of the subtracted spectral density located at m S2 = 1.50 GeV points to the f 0 (1500) meson [29]. Since f 0 (1500) has a narrower width 112 MeV [12], it cannot accommodate the broad width shown in Fig. 9 alone. Hence, it is likely that f 0 (1370) and f 0 (1710) also have gluonic contents and contribute to the spectral density, in accordance with the prevailing consensus in the literature [33]. That is, all f 0 (1370), f 0 (1500) and f 0 (1710) are glue-rich states (nevertheless, a recent phenomenological analysis on J/ψ radiative decays prefers a higher scalar glueball mass resulting from the mixing with heavier scalar states [95,96]). The shorter peak located at m S1 ≈ 0.60 GeV might arise from the contribution of the f 0 (500) meson with a broad width. The smaller area under the peak implies a lower f 0 (500) decay constant defined via the gluon field, and less gluonic content in this light scalar, consistent with the observation in [31]. This implication does not depend on the quark structure of the f 0 (500) meson [97]. Note that a little amount of gluonium in the f 0 (980) meson with a narrow width 10 MeV cannot be excluded, as hinted by the spectral densities in Fig. 8. The above interpretation agrees with the conclusion drawn based on a five-state-mixing scenario in a nonlinear chiral Lagrangian framework [32,98] and with the analysis in [13]. The mass gaps between f 0 (500) and f 0 (980), and among f 0 (1370), f 0 (1500) and f 0 (1710) are too small to be resolved by our method with the characteristic scale Λ ∼ O(1) GeV 2 , so only two peaks are revealed in the spectral density.

C. Pseudoscalar Glueball Masses
Next we extract the pseudoscalar glueball mass by solving Eq. (38) from the corresponding set of OPE inputs in Eq. (33) with the triple-gluon condensate in Eq. (40). The prescription for fixing the number N of polynomials in an expansion is the same, and the plots for ρ P (s) are similar to those in Fig. 8. The behavior of the subtracted spectral density ∆ρ P (s, Λ) for Λ = 1.5 GeV 2 with the expansion up to N = 14 generalized Laguerre polynomials L (2) n (y) is exhibited in Fig. 11, from which we read the pseudoscalar masses. The coefficients a n in this case have not yet grown quickly, but the positivity constraint forces us to terminate the polynomial expansion at N = 14. The maximal number N increases with Λ under the positivity constraint till Λ = 2.3 GeV 2 , above which the spectral density is always positive, but the matrix elements of M −1 and the coefficients a n go out of control as N > 22. Therefore, we stick to N = 22 for the polynomial expansion in the range Λ = 2.3-3.5 GeV 2 . Likewise, we observe the clear double-peak structure in Fig. 11 with the locations m 2 P1 ≈ 0.5 GeV 2 and m 2 P2 = 3.1 GeV 2 , which are above the scalar ones m 2 S1 ≈ 0.4 GeV 2 and m 2 S2 ≈ 2.3 GeV 2 in Fig. 9, respectively. We scan the pseudoscalar masses m P1 and m P2 in the range 1 GeV 2 < Λ < 3.5 GeV 2 , and display their dependencies on Λ in Fig. 12, where the two curves describe the two peak locations of the subtracted spectral density ∆ρ P (s, Λ). It is found that the lower and upper curves almost remain flat in the interval 1.0 GeV 2 < Λ < 2.3 GeV 2 , around m P1 = 0.71 GeV and m P2 = 1.75 GeV, and then go up monotonically as expected. The stability window is wider than in the scalar glueball case illustrated in Fig. 10. We estimate the theoretical errors in our method using the extreme values in the window, and get m P1 = (0.71 ± 0.02) GeV and m P2 = (1.75 ± 0.02) GeV, whose tiny errors reflect the remarkable stability of our solutions. We investigate the theoretical uncertainties from the variation of the involved parameters in a similar manner. The decrease of the strong coupling α s from 0.5 to 0.4 yields only 7% enhancement on the determined pseudoscalar masses, which justifies the neglect of the renormalzation-group evolution in our analysis. The typical ±20% change of the gluon condensate α s G 2 causes about ±5% effect. The choice for That is, the theoretical uncertainties are under control in our formalism. The major peak of the subtracted spectral density located at m P2 ≈ 1.75 GeV may point to the η(1760) meson, whose broad width about 240 MeV [12] can accommodate the width in Fig. 12 alone. This mass is just a bit higher than the scalar glueball one determined in the previous subsection, as anticipated from the minor difference between their OPE inputs and argued in [99]. Our result is lower than most predictions in the literature, which are above 2 GeV, as stated in the introduction. We stress again that measurements of J/ψ radiative decays do not confirm any glue-rich pseudoscalar resonances with masses greater than 2 GeV [12]. The branching ratio of the J/ψ → γX(2370) decay is so small [12], that X(2370) is unlikely to be the pseudoscalar glueball, even if it carried the correct quantum numbers. On the contrary, the η(1760) meson is abundantly produced in J/ψ radiative decays, but not seen in the J/ψ → γγV decays [50]. Another nearby pseudoscalar X(1835) is produced less abundantly in J/ψ radiative decays, and seen in the J/ψ → γγφ decay [101]. It should be reminded that the relevant experimental studies are not yet conclusive. For instance, it is puzzling that the J/ψ → γ(η(1760) →)ωω branching ratio, despite of the stronger phase space suppression, is about one order of magnitude larger than the J/ψ → γ(η(1760) →)ρ 0 ρ 0 one [12]. Nevertheless, the currently available data do support that η(1760) is a glueball.
Our solution disfavors the speculation, deduced from a pseudoscalar meson mixing formalism based on the anomalous Ward identity [51][52][53], that the η(1405) is the lightest pseudoscalar glubeball [102]. It should be pointed out, however, that a pseudoscalar glueball mass as heavy as 1.75 GeV is not excluded in [51], when some inputs are allowed to vary. A similar mixing formalism with more conservative assumption [55] shows that the pseudoscalar glueball tends to be heavier. Our prediction 1.75 GeV matches their results with a large angle φ G for the mixing between the pure glueball and the flavor-singlet light quark states. It is fair to allege, based on the theoretical uncertainties, that the η(1405/1475) mesons [119] contain some gluonium components, which, though, are not dominant as found in [13]. Indeed, the η(1405/1475) mesons are produced copiously in J/ψ radiative decays, and also seen in the J/ψ → γγρ decay [12].
The shorter peak located at m P1 ≈ 0.71 GeV between the η and η meson masses comes from the combined contributions of these two states with similar weights. A low-lying state with mass around 1 GeV has been also identified and assigned to the η meson in the lattice calculation [100], when the topological charge density with a strong coupling to flavor-singlet light quark states, the same as in the present work, is employed to define the correlation function. This observation is in accordance with the implementation of the η-η -glueball mixing, which has been intensively discussed in the literature. The comparison of the relative heights between the two peaks in Figs. 9 and 11 suggests that η and η have more gluonic content than f 0 (500) and f 0 (980) do. Our solutions hint that a more complete mixing scenario involving η, η , η(1405/1475) and η(1760) is needed for a deeper understanding of the pseudoscalar glueball properties. We will examine whether more quantitative information on the scalar and pseudoscalar mixings can be drawn from our formalism in the future. For the same reason, the mass gap between η and η is too small to be resolved by our method with the characteristic scale Λ ∼ O(1) GeV 2 , although they have quite narrow widths.

D. Tensor Glueball Mass
We extend the above formalism to the study of the tensor glueball (2 ++ ), for which the correlation function is defined as with η µν = g µν − q µ q ν /q 2 , and being the energy-momentum stress tensor of QCD. The OPE of the correlation function Π T (q 2 ) in the deep Euclidean region of q 2 is given by [14] Π OPE up to the dimension-eight condensate, ie., to the power correction of 1/(q 2 ) 2 , which can be approximated under the vacuum factorization assumption by Repeating the same steps, we notice that the tensor glueball case differs much from those of the scalar and pseudoscalar glueballs. The positivity requirement on the spectral density forces us to terminate the expansion in terms of the generalized Laguerre polynomials L (2) n (y) at very low N : the maximal N 's are found to be 3, 4, 6 and 8 for Λ = 1.5, 2.0, 3.0 and 4.0 GeV 2 , respectively. The maximal number N for a polynomial expansion increases with Λ under the positivity constraint, but it reaches only N = 10 even when Λ is as high as 5.0 GeV 2 . Hence, we accept the expansions up to fewer polynomials in the analysis of the tensor glueball mass, and obtain the spectral density ρ T (s) and the subtracted spectral density ∆ρ T (s, Λ) for Λ = 3.0 GeV 2 and N = 6 given in Fig. 13. The coefficients a n in the expansion have not increased rapidly, with a 6 /a 5 ≈ 1.2. If N = 7 is chosen, ρ T (s) will become negative and violate the positivity constraint in the range 0 < s < 0.7 GeV 2 . No evident bump like those in Fig. 7 appears in the curve for ρ T (s), which always ascends in the range 1.0 GeV 2 < Λ < 5.0 GeV 2 . After the huge perturbative background is removed, the subtracted spectral density ∆ρ T (s, Λ) exhibits a single peak with diminishing height (note the scale on the vertical axis in Fig. 13). The single-peak structure of ∆ρ T (s, Λ) persists till Λ = 5.0 GeV 2 , and becomes broader with Λ. It has been underlined that a peak of the spectral density cannot be interpreted as a physical state, unless a stability window in Λ for its location exists. We thus investigate the dependence of the peak location on Λ in the range 1.5 GeV 2 < Λ < 5.0 GeV 2 , and present it in Fig. 14. It turns out that the mass m T always grows with Λ monotonically, namely, a stable region for the tensor glubeball mass is not found. The number of the generalized Laguerre polynomials is smaller than 3 for Λ < 1.5 GeV 2 , but the resultant masses, if depicted in Fig. 14, follow the same trend. The curve is a bit bumpy, because the polynomials are few, such that the discontinuity at each increment of N is more significant. It has been elaborated [2] that at least two condensates with different dimensions are necessary for establishing the ρ meson state. Hence, we speculate that the absence of a resonance solution for the tensor glueball may be ascribed to the insufficient nonperturbative input in the present setup: the single four-gluon condensate may not be able to fix the tensor glueball mass. An OPE, which includes higher-dimensional operators compared to Eq. (47), is urged. Note that a minimum of m T ≈ 2.0 GeV was extracted from the ratio of two selected moments by varying the Borel mass in sum rules [13]. However, no minimum existed as the the continuum threshold was varied, a situation different from the scalar and pseudoscalar glueball cases [13]. We suspect that more moments considered in our formalism may have imposed a stronger constraint on the existence of the tensor glueball. This issue is worth more thorough studies.

IV. CONCLUSION
We have refined our previous proposal for handing QCD sum rules as an inverse problem by solving dispersion relations with OPE inputs directly. A nonperturbative spectral density is expanded in terms of a suitable set of generalized Laguerre polynomials up to some degree, according to its boundary condition at vanishing energy. A dispersive integral and power corrections to an OPE from condensates render possible expansions of both sides of a dispersion relation into power series in 1/q 2 . The coefficients in the power series on the hadron side can then be derived in the inverse matrix method from the known OPE coefficients on the quark side straightforwardly. An additional polynomial in the expansion for the spectral density appears as a higher-power correction in 1/q 2 to the dispersion relation, such that the convergence of the OPE guarantees the convergence of the polynomial expansion for the spectral density. Certainly, a solution from an ill-posed inverse problem will go out of control, as the matrix dimension becomes sufficiently large. However, a reasonable solution can be obtained in most cases, as elaborated in Sec. IIA by means of mock data constructed from several test functions. The employment of the generalized Laguerre polynomials introduces an arbitrary scale Λ into the theoretical framework, which characterizes the resolution power of the method. It has been explained that a peak in the spectral density can be identified as a physical resonance, if its location is stable against the variation of Λ.
The applications of the above formalism to the determinations of the ρ meson, scalar glueball and pseudoscalar glueball masses have been quite successful in the sense that the stability windows in Λ do exist. By producing the ρ meson mass m ρ = 0.78 GeV, we fixed the gluon condensate α s G 2 and the factorization violation parameter κ, which are both within the ranges commonly accepted in the literature. The double-peak structures, appearing naturally in the spectral densities for the scalar and pseudoscalar glueballs, indicate some gluonium components in the light quark states. The locations of the major broad peaks, 1.50 and 1.75 GeV, point to the f 0 (1500) and η(1760) as the glue-rich states in the scalar and pseudoscalar sectors, respectively. The former, together with f 0 (1370) and f 0 (1710), have been well recognized and extensively analyzed in mixing frameworks. However, the latter, with its mass below most predictions from quenched lattice QCD and sum rules, has not, and deserves more theoretical and experimental endeavors. We stress that we did not find a resonance solution for the spectral density associated with the tensor glueball. As shown in our previous work, at least two condensates with different dimensions in the OPE input are required for establishing the ρ meson state. We have speculated that the absence of a solution for the tensor glueball may be due to the insufficient nonperturbative input at present. An OPE for the corresponding correlation function with higher-dimensional operators is thus in demand.
A merit of our approach is that the correlation functions for the scalar and pseudoscalar glueballs at zero momentum can be calculated from the subtracted spectral densities introduced in this paper, which realize explicitly the ultraviolet regularization required for dispersive integrals in the literature. The former served to discriminate the lattice estimate for the triple-gluon condensate from the others, and the latter gave rise to our prediction for the topological susceptability χ 1/4 t = 75-78 MeV, whose range overlaps with that from lattice QCD. Combined with the other findings on the measured ρ meson mass and the widely accepted scalar glueball mass around 1.5 GeV in the same formalism, and the experimental observations from J/ψ radiative decays, we tend to advocate that the η(1760) meson is a promising candidate for the pseudoscalar glueball.
We have explored various sources of theoretical uncertainties in our method, which are all under control. The error from the variation of Λ in the stability window is tiny, as reflected by the flat curves for the resonance masses. The error from the renormalization-group evolution is also minor, because it is the relative importance among the different terms in an OPE that determines the resonance masses, which is not sensitive to the running effect. The variations of the condensates are not crucial either: a typical 20% change of the gluon condensate causes about 5% effects on the ρ meson and glueball masses. Compared to conventional sum rules, outcomes from directly solving dispersive relations have less model dependence. Besides, our theoretical framework can be improved systematically by taking into account higher-order and higher-power corrections to an OPE on the quark side. Note that the dimension-eight condensates are still quite uncertain [106][107][108][109], on which more progress is needed.
Our work suggests that the multiple-resonance mixing scenario for the pseudoscaalr mesons should include η(1760) for completeness. To explore mixing properties in our formalism, one has to consider additional off-diagonal correlation functions, in which one of the gluon operator is replaced by a quark one [110,111]. The decay constants of glueballs ought to be extracted in mixing frameworks, since their definitions depend on the currents adopted. Moreover, the input from the triple-gluon condensate, which affects mixing patterns [36], should be fixed accurately first. Though our method is powerful for studying properties of low-lying resonances, it is difficult to probe excited states and finer structures, such as the ρ-ω mixing. To attempt the former, one may resort to the multiple-pole parametrization plus an arbitrary continuum contribution for a spectral density as in conventional sum rules [94,[112][113][114][115][116]. However, it is likely that we can avoid the ad hoc prescriptions for choosing relevant hadronic parameters, such as continuum thresholds [117]. It is worthwhile to extend our formalism to investigations of three-gluon states [118] and other hadronic states. There is no doubt that that our proposal will have wide applications to analyses of nonperturbative observables.