Mass Calculations of Light Quarkonium, Exotic $J^{PC}=0^{+-}$ Hybrid Mesons from Gaussian Sum-Rules

We extend previous calculations of leading-order correlation functions of spin-0 and spin-1 light quarkonium hybrids to include QCD condensates of dimensions five and six, with a view to improving the stability of QCD sum-rules analyses in previously unstable channels. Based on these calculations, prior analyses in the literature, and its phenomenological importance, we identify the exotic $J^{PC}=0^{+-}$ channel as the most promising for detailed study. Using Gaussian sum-rules constrained by the H\"older inequality, we calculate masses of light (nonstrange and strange) quarkonium hybrid mesons with $J^{PC}=0^{+-}$. We consider single narrow, single wide, and double narrow resonance models, and find that the double narrow resonance model yields the best agreement between QCD and phenomenology. In both non-strange and strange cases, we find hybrid masses of $2.60$ GeV and $3.57$ GeV.


Introduction
It has long been conjectured that hadrons could exist beyond the conventional quark model of quark-antiquark (qq) mesons and three-quark (qqq,qqq) baryons. In particular, colour-singlet hybrid mesons consisting of a quark, antiquark, and explicit gluonic degree of freedom have a long history [1]. While evidence of hadronic structures outside of the conventional model has been accumulating with experimental observations and confirmations of tetraquarks [2,3,4] and pentaquarks [5], an experimental confirmation of hybrid mesons has eluded observation. Designed to search for light hybrid mesons (particularly those with exotic J P C that do not exist in the conventional quark model), the GlueX experiment at Jefferson Lab [6] is currently underway, and is anticipated to give crucial insight into the existence and structure of light hybrids.
In [33,34,35,36], it was found that the inclusion of higher-dimension condensates stabilized previously unstable LSRs analyses from [11] of hybrids containing heavy quarks. Therefore, we provide a systematic computation of leading-order (LO) 5d and 6d condensate contributions for all light quarkonium hybrids of spin-zero and spin-one. Unfortunately, these higher-dimension condensates do not stabilize the unstable light hybrid LSRs analyses as they do for heavy hybrids. However, in [11], it was proposed that the instability in the LSRs might be resolved by accounting for finite width effects, an issue also raised in [13]. To explore width effects and the possibility of excited states, we depart from previous LSRs methods. Gaussian sum-rules (GSRs) [37] are sensitive probes of width effects and both ground and excited states, and have been shown to be a powerful and versatile analysis methodology [38,39,40,41]. Thus, in this article, we use GSRs to investigate the possibility of distributed resonance strength in the exotic 0 +− light hybrid channel.
In Section 2, we calculate LO spin-0 and spin-1 correlation functions of light quarkonium hybrid currents, including condensates up to 6d. We discuss the higher-dimension condensates' lack of effect on the stabilization of LSRs analyses. We conclude that the most promising avenue for further exploration is a GSRs analysis of the exotic 0 +− channel. Section 3 reviews the GSRs formalism, and a theoretical constraint on the GSRs based on the Hölder inequality is developed in Section 4. The GSRs analysis methodology and results for the 0 +− channel are presented in Section 5 with concluding remarks in Section 6.

Hybrid Currents and Correlation Functions
To investigate light quarkonium hybrids, we use currents of the form where q is a light (nonstrange or strange) quark field and t a are generators of the fundamental representation of SU (3). Each combination of G a µν ∈ G a µν ,G a µν = 1 2 µνρσ G a ρσ and Dirac structure Γ ν together corresponds to particular values of J P C [11,15]; these combinations are summarized in Table 1.  (1).
For each current (1), we calculate and decompose a diagonal correlation function as follows: where Π (0) probes spin-0 states and Π (1) probes spin-1 states. The calculation of (2) is performed in the framework of the operator product expansion (OPE), In (4), the vacuum expectation value (VEV) of a time-ordered, non-local product of composite operators is expanded in a series, each term of which is a product of a perturbative Wilson coefficient C n (x) and a nonzero VEV of a local composite operator O n (0), i.e., a condensate. The condensates parameterize the nonperturbative nature of the QCD vacuum, and we include in our correlator calculations the following set: respectively the 3d quark condensate, the 4d gluon condensate, the 5d mixed condensate, the 6d gluon condensate, and the 6d quark condensate. In (5)-(9), superscripts on quark fields are colour indices whereas subscripts are Dirac indices and σ µν = i 2 [γ µ , γ ν ]. Regarding Wilson coefficients, we consider LO calculations in α s , and we compute O m 2 light quark mass corrections to perturbation theory as a way to distinguish between the nonstrange-and strange-flavored cases, similar to [42]. Also, the values of (5), (7), and (9) depend on whether the light quarks are nonstrange or strange. The diagrams representative of the correlation function calculation are displayed in Figure 1. We use dimensional regularization in D = 4+2 dimensions at MS renormalization scale µ. The program TARCER [44] is utilized to reduce the (a) Diagram I (LO perturbation theory) Figure 1: The Feynman diagrams calculated for the correlator (2). Feynman diagrams were created using JaxoDraw [43]. Table 2: Coefficients of the logarithmic and divergent terms of the perturbative and condensate contributions to the correlation function (10) for the J P C summarized in Table 1.
resulting integrals to a selection of well-known master integrals using the Tarasov recurrence relations [45,46]. All of the correlators defined between (1) and Table 1 can be written in general as where we have suppressed the superscript (J) on the left-hand side. The coefficients A i and B j contained in (10) are given in Tables 2 and 3 respectively. We note that, as Diagram IV has no loops, A 5 is trivially zero. In all channels, perturbation theory, the 3d quark condensate term, and the 4d gluon condensate term were benchmarked against [11]. The 0 −− and 1 −+ correlators were benchmarked against [15].

QCD Sum-Rules
Each function Π (J) (q 2 ) defined in (3) satisfies a dispersion relation at Euclidean momentum where we have again suppressed the superscript (J). In (11), t 0 is a hadron production threshold and · · · are subtraction constants, together a third degree polynomial in Q 2 . Equation (11) connects theoretical predictions of QCD, i.e., Π(Q 2 ) on the left-hand side, to properties of hadrons contained in ImΠ(t), the hadronic spectral function, on the right-hand side. Table 3: Coefficients of the finite terms of the perturbative and condensate contributions to the correlation function (10) for the J P C summarized in Table 1. Regarding (11), to eliminate subtraction constants and to accentuate the low-energy region of the integral on the right-hand side, some transformation is typically applied. A popular choice is to formulate unsubtracted LSRs of (usually nonnegative) integer weight k, at Borel parameter τ [47,48,49,50]. Details on how to evaluate (12) for a correlator such as (10), denoted Π QCD from here on to emphasize that it is a quantity calculated using QCD, can be found in the literature (e.g., [47]). The result is for k ∈ {0, 1, 2, . . .} and where Recall, the A i are given in Table 2.
In (11), we impose on ImΠ(t) a general resonances-plus-continuum model with onset of the QCD continuum at threshold s 0 , where ρ had (t) represents the resonance content of the hadronic spectral function and θ(t) is the Heaviside step function. To isolate the resonance contributions to the LSRs, we consider (continuum-) subtracted LSRs Then, Equations (11)-(13), (15), and (16) together imply that and (again) ImΠ QCD (t) is given in (14).
There are a number of interesting observations we can make concerning the LSRs of light quarkonium hybrids. In particular, the 6d gluon condensate terms do not contain a logarithm, i.e., A 6 = 0 for all J P C values considered (see Table 2), and hence do not contribute to the imaginary part (14). This result is surprising: both Diagrams V and VI (see Figure 1) have logarithmic contributions, but they cancel when the two diagrams are added together. Thus, the LO 6d gluon condensate terms cannot stabilize light quarkonium hybrid LSRs analyses as they have done in some heavy quarkonium hybrid analyses [33,34,35,36].
Another observation relates to the mixed condensate contributions. Using (12), if we try to formulate k = −1 (i.e., lower-weight) unsubtracted LSRs, we get a piece that formally looks like the right-hand side of (13) at k = −1 and another piece: If A 7 = 0, then neither piece is well-defined: the integral from (13) diverges and (19) contains a −1 field theory divergence. But for J P C ∈ {1 −+ , 1 ++ }, we find that A 7 = 0 which allows for the construction of lower-weight LSRs in these two channels. Unlike the k = 0 LSRs, the k = −1 LSRs do receive contributions from the 6d quark and gluon condensates as both B 5 and B 6 are nonzero. An analysis of these k = −1 LSRs does require some knowledge of the subtraction constants in (11). As noted in Section 1, in the multi-channel LSRs analysis of [11], the 0 +− , 0 −+ , 1 +− , and 1 −+ sectors were unstable. The 1 −+ has since been stabilized using lower-weight LSRs [15], and the 0 −+ has been stabilized [16] using a different current than that used in [11]. That leaves the non-exotic 1 +− and the exotic 0 +− channels. Given the GlueX emphasis on exotics and the possible complicated features of mixing between hybrids and conventional quark mesons in the 1 +− channel, we focus our attention on 0 +− light quarkonium hybrids. Attempts to stabilize the 0 +− channel have involved higher-dimension currents combined with lowerweight sum-rules requiring estimation of the dispersion-relation low-energy constant within the analysis [13]. Because higher-dimension currents tend to enhance the continuum, the mass determination combined with an estimated low-energy term merits further study.
As in [11], we perform a conventional single narrow resonance (SNR) LSRs analysis of the 0 +− channel by letting in (17) where f is the resonance coupling and m H is its mass. We include our higher-dimension condensate contributions as well as updated QCD parameter values, yet the analysis remains unstable. The 5d mixed condensate term in the LSRs is small, and, as noted above, the 6d condensates do not contribute at all. In [11], it was suggested that the instability in this channel could be related to a distribution of resonance strength. To investigate this possibility, we use GSRs, an alternative to LSRs which provide a fundamentally different weighting of the hadronic spectral function that makes them well-suited to analyzing distributed resonance strength hadron models. Unsubtracted GSRs of integer weight k are defined as [37] G Details on how to evaluate (21) for a correlator such as (10) can be found in [37,38,39]. The result is for k ∈ {0, 1, 2, . . .} and where 1 π ImΠ QCD (t) is given in (14). Subtracted GSRs are defined in much the same way as subtracted LSRs leading to the following GSRs analogues of (17) and (18): where The difference between (17)- (18) and (23)- (24) is, of course, the kernel of the integrals: a decaying exponential for LSRs and a Gaussian for GSRs. From (23), we see that GSRs smear ρ had (t) over a Gaussian distribution of width √ 2τ centered atŝ. In the τ → 0 + limit, we have which, when applied to (23), yields Hence, at least in principle, ρ had (t) can be extracted directly from GSRs. Realistically, the τ → 0 + limit cannot be achieved, however, because, through renormalization-group (RG) improvement (see Section 5), the renormalization scale at which α s is evaluated decreases with decreasing τ [37]. Nevertheless, it is desirable to use low values of τ to minimize the smearing of ρ had (t) by the kernel of the GSRs. To further emphasize this, we draw upon an analogy introduced in the seminal GSRs paper [37]. Gaussian sum-rules satisfy the classical heat equation reinterpreting the parameterŝ as "position", the Gaussian width τ as "time", and the GSRs G k (ŝ, τ, s 0 ) as "temperature". The smaller the value of τ (i.e. the less "time" that has passed), the better we may assess the original (i.e., τ → 0 + ) "temperature" distribution (i.e., s k 1 π ρ had (ŝ)). Compared to LSRs, GSRs permit greater access to the structure of ρ had (t). The LSRs methodology is specifically formulated to accentuate the ground state while suppressing excited states. With GSRs, this need not be the case asŝ, the position of the Gaussian kernel's peak, is a free parameter. By varyingŝ, GSRs can probe excited states with the same sensitivity as the ground state. As such, GSRs are generally preferable to LSRs when studying distributed resonance strength models.
Integrating (23) with respect toŝ gives from which we recognize the quantity on the left-hand side as the finite-energy sum-rule of weight k. As shown in [37], a resonance plus continuum model evolved through the diffusion equation only reproduces the QCD prediction at large energy scales if s 0 is constrained by (28).

Hölder Inequality
Previous investigations of hadronic systems using LSRs have employed Hölder inequalities to restrict the set of allowed τ and s 0 values [51,52,53]. The Hölder Inequality is expressed generally as 1 q under the condition 1 p and where dµ is an arbitrary integration measure. From positivity of the hadronic spectral function for diagonal correlators, we can use ImΠ QCD (t) > 0 to form an integration measure. Substituting this integration measure into (32) leads to restrictions on the allowed values of s, τ , and s 0 in the GSRs. We consider the inequality (32) with the assignments dµ = ImΠ QCD (t)dt (34) where α + β is a non-negative integer. Defining the inequality (32) becomes where we have used G k (τ,ŝ, s 0 ) > 0, the weakest constraint on the GSRs that emerges from positivity of the spectral function. We define ω as follows: and consider (40) with zero-weight GSRs (i.e., α = β = 0), Equations (38), (39), and (41) together imply that which, when substituted into (42), gives Following [53], we set which implies We can perform a local analysis of (47) by expanding about δτ = 0, where primes indicate τ -derivatives. Then, (41) and (48) together imply (49) At some (τ,ŝ, s 0 ), if the GSR G 0 (ŝ, τ, s 0 ) is to be consistent with a positive hadronic spectral function, then it must satisfy the inequality (47).

Analysis Methodology and Results
Before we can analyze 0 +− light quarkonium hybrids using (31), we need to discuss the QCD parameters appearing in (10), i.e., the coupling, the quark mass, and the QCD condensates.
To implement RG improvement we replace α s and m in (10) by one-loop, MS running quantities [37]. In our analysis, we use the QCD running coupling anchored at the τ -lepton mass, where we use PDG [54] values for the τ mass and For the light quark masses, we use where for nonstrange quarks and m(2 GeV) = 96 +8 −4 MeV for strange quarks [54]. In both (50) and (52), we set n f = 4.
Renormalization-group arguments identify our renormalization scale as µ = τ 1/4 [37,38], putting a lower bound on our choice of τ restricted by the reliability of perturbation theory. We therefore restrict our analysis to τ ≥ M τ , approximately equivalent to τ = 10 GeV 4 . We work with an upper bound τ ≤ 20 GeV 4 , as discussed in detail in Section 5.
Turning to the condensates, the value of the RG-invariant quantity mqq is well-known from PCAC [55]. Using the conventions of [56], we have where PDG values are used for the meson masses [54] and the decay constants are [57] We use the following value for the 4d gluon condensate [58]: The nonstrange-and strange-flavored 5d mixed condensates are estimated by [59,60] to be Finally, we note that while the 6d quark and gluon condensates were included in the correlator calculation (10), Table 2 shows that neither contributes to the k = 0 GSRs (24) or NGSRs (29). As noted in Section 3, a SNR analysis of 0 +− light quarkonium hybrids fails, and so we turn our attention to models with distributed resonance strength. First, we consider the quantity [39] where the moments, M k,n (τ, s 0 ), were defined in (30). Combining (23) and (59) gives For a SNR model, substituting (20) into (60) yields In Figure 2, we plot σ 2 0 (τ, s 0 ) − 2τ versus τ for nonstrange quarks at several values of s 0 over the range 10 GeV 2 ≤ s 0 ≤ 30 GeV 2 . Clearly, σ 2 0 (τ, s 0 ) − 2τ = 0, in violation of (61), which provides further motivation to consider models other than the SNR. An analogous analysis for strange quarks leads to the same conclusion.
If the distributed resonance strength indicated by Figure 2 is due to a single wide resonance (SWR), then we can determine a rough lower bound on the resonance's width using a rectangular pulse resonance model, where f is the resonance's coupling, Γ is its width, and m H is its mass. Substituting (62) into (60) gives From (64), we see that Γ decreases as m H increases. However, to ensure that the resonance does not merge with the continuum, we require which implies that the largest possible resonance mass for a particular s 0 is given by where we have used (63). By letting m H → m H, max in (64), we find that the smallest possible resonance width for a particular s 0 is given by From Figure 2, we see that σ 2 0 (τ, s 0 ) − 2τ shows almost no τ -dependence; hence, the same can be said about Γ min (τ, s 0 ). In Figure 3, we plot Γ min (τ, s 0 ) versus s 0 at τ = 10 GeV 4 for nonstrange quarks. An analogous plot for strange quarks looks nearly identical. At s 0 = 10 GeV 2 , we find that Γ min ≈ 1.46 GeV, far larger than a typical hadron width. As s 0 increases, so too does Γ min . For these reasons, we abandon SWR models in favour of a multi-resonance model. We consider a double narrow resonance (DNR) model where f 1 , f 2 and m 1 , m 2 are the resonances' couplings and masses respectively. Substituting (68) into (31) gives At fixed values of τ and s 0 , we perform a fit of (69) overŝ 1 to find best fit parameters for r, m 1 , and m 2 . In Figure 4, we plot the best fit r versus s 0 at τ = 10 GeV 4 for nonstrange quarks. Again, an analogous plot for strange quarks looks nearly identical. From the s 0stability in r versus s 0 , we determine an optimized continuum onset for both the nonstrangeand strange-flavored cases as s opt 0 = (14.5 ± 1.2) GeV 2 where the uncertainties originate from the QCD input parameters; details of the uncertainty analysis are discussed below. Then, a fit to (69)  in the strange-flavored case. Figure 6 shows comparisons between the the NGSRs and the DNR model (respectively the left-and right-hand sides of (69)) for parameters (72)-(74) at τ = 10 GeV 4 and τ = 20 GeV 4 . We note that the strange and nonstrange 0 +− hybrid mass predictions are degenerate within the uncertainties of our analysis. Utilizing the Hölder Inequality test (49), we can perform a consistency check on our analysis. To determine whether (49) is satisfied within the expected uncertainties of the GSRs, we examine the inequality at s opt 0 = 14.5 GeV 2 for various values of τ . Because our QCD calculations of Wilson coefficients are truncated perturbative series in α s , in addition to the QCD parameter uncertainties, we use the 1 −+ channel [61] to provide an estimated next-order  perturbative correction characteristic of hybrid correlators. We find that the Hölder inequality constraint (49) is violated for τ 20 GeV 4 , and the inequality test for the minimum value τ = 10 GeV 4 is shown in Figure 7. Thus, the τ range used in our analysis, 10 GeV 4 < τ < 20 GeV 4 , is consistent with the Hölder Inequality. We verify the s 0 optimization (71) obtained from Figure 4 by looking at an independent analysis developed in [39] based on the properties of theŝ peak position (maximum) of the NGSRs. For the SNR model (20) theŝ-peak occurs atŝ = m 2 , independent of τ . Thus, an alternative s 0 -optimization criterion for the SNR is minimizing the τ -dependence of the peak positionŝ peak (τ, s 0 ) defined implicitly from For the DNR model (68), the peak position acquires τ -dependence modeled bŷ   By minimizing (80) with respect to A, B, C, D, and s 0 , we find an optimum continuum threshold s opt 0 = 14.0 GeV 2 in excellent agreement with the value obtained in (71). To obtain errors in s opt 0 , r, m 1 , and m 2 , we examine how the errors in the QCD parameters impact the values of these optimized parameters by varying each independently and examining the impact on the model parameters. Additionally, there exists a methodological error in determining s opt 0 as the variance in the QCD parameters will affect the stability point of r.
Contributions to the error in s opt 0 are summarized in Table 4 and contributions to the error in the DNR model parameters are summarized in Tables 5-7. The dominant error in s opt 0 comes from the variation in αG 2 ; in determining errors in the DNR parameters, the error in r is driven by the variation in αG 2 while the dominant errors in the masses m 1 and m 2 arise from variations in s opt 0 , followed by αG 2 . Errors in qq and gqσGq contribute negligibly in the error of all DNR parameters. Adding the values summarized in Tables 4-7 in quadrature gives us a conservative error estimate summarized in Table 8; as the driving errors in each parameter are approximately equivalent for the upper and lower bounds of the corresponding QCD parameters, we express our DNR parameters (71)-(77) with symmetric error, taking the most conservative bound.

Discussion
We have calculated 5d and 6d QCD condensate contributions to all spin-0 and spin-1 light quarkonium hybrid correlators with the goal of obtaining QCD LSRs mass predictions in the previously-unstable channels of [11]. However, the 6d gluon and quark condensate contributions do not have an imaginary part and hence do not contribute to the LSRs. Also, the 5d mixed condensate contributions turn out to be small. We therefore focused on the suggestion of References [11,13] that a distribution of resonance strength could be the source of instability, a scenario ideally suited to GSRs methods [37,38,39]. The 0 +− channel was chosen for detailed investigation because of its phenomenological significance in light of the GlueX experiment.
In examining the SWR (62) and DNR (68) models, we found that the DNR model provided excellent agreement between QCD and phenomenology. (See Figure 6.) The SWR model was rejected on the basis of an atypically large resonance width. In the DNR model, we find degenerate predictions in the case of both nonstrange and strange quark states from the 0 +− current: a 2.60 ± 0.14 GeV state (2.60 ± 0.14 GeV in the strange case) with 29% relative coupling, and a state at 3.57 ± 0.15 GeV (3.57 ± 0.13 GeV) with 71% relative coupling. The smaller coupling of the light state suggests the possibility of mixing with a tetraquark because the expected tetraquark mass range is above 2 GeV [62].
The lighter state is consistent with recent lattice results that find a predominantly nonstrange state around 2.4 GeV and a predominantly strange state around 2.5 GeV in the 0 +− channel with m π = 391 MeV [25]. Our lighter-state mass determination is somewhat larger than the 2.1-2.5 GeV range of central values in [13]. The literature does not provide a clear interpretation of the heavier 0 +− state; however lattice studies [25] of the lightest hybrid meson supermultiplet suggest that the 0 +− state exists as part of an excited hybrid supermultiplet with radially-excited qq pair (i.e., quark total angular momentum L qq = 1). We suggest that this heavier second state arising in the GSRs is a manifestation of an excited hybrid state.
In conclusion, we investigated light quarkonium, exotic (J P C = 0 +− ) hybrid mesons with SWR and DNR models using a GSRs analysis. We disfavoured the SWR model as the predicted resonance width was too large. The double-narrow resonance model yielded two 0 +− hybrid states: (2.60±0.14) GeV and (3.57±0.15) GeV ((2.60±0.14) GeV and (3.57±0.13) GeV in the strange case). Additionally, we explored using the Hölder inequality derived for the GSRs as a consistency check on our analysis.