$1^{-+}$ Hybrid in $J/\psi$ Radiative Decays from Lattice QCD

We present the first theoretical prediction of the production rate of $1^{-+}$ light hybrid meson $\eta_1$ in $J/\psi$ radiative decays. In the $N_f=2$ lattice QCD formalism with the pion mass $m_\pi\approx 350$ MeV, the related electromagnetic multipole form factors are extracted from the three-point functions that involve necessarily quark annihilation diagrams, which are calculated through the distillation method. The partial width of $J/\psi\to \gamma \eta_1$ is determined to be $2.29(77)~\mathrm{eV}$ at the $\eta_1$ mass $m_{\eta_1}=2.23(4)$ GeV. If $\eta_1$ corresponds to the recently observed $\eta_1(1855)$ in the process $J/\psi\to \gamma\eta_1(1855)\to \gamma \eta\eta'$ by BESIII, then the branching fraction $\mathrm{Br}(J/\psi\to \gamma\eta_1(1855))$ is estimated to be $6.2(2.2)\times 10^{-5}$, which implies $\mathrm{Br}(\eta_1(1855)\to \eta\eta')\sim 4.3\%$.

We present the first theoretical prediction of the partial decay width of the process J/ψ → γη1, where η1 is the lightest flavor singlet 1 −+ hybrid meson. Our N f = 2 lattice QCD calculation at mπ ≈ 350 MeV results in the η1 mass mη 1 = 2.23 (4) GeV and the related electromagnetic form factors M1(0) = −4.73(74) MeV, E2(0) = 1.18 (22) MeV, which give Γ(J/ψ → γη1) = 2.04(61) eV. These form factors can be applied to the physical N f = 3 case where there should be two hybird mass eigenstates η (l) 1 and η (h) 1 due to the singlet-octet mixing. It is shown that the ratio of the branching fractions Br(J/ψ → γη (l,h) 1 → γηη ) is inversely proportional to the ratio of the total widths of η (l,h) 1 . Given our results and the mixing angle derived by a previous lattice study, whether η1(1855) is assigned to be η

I. INTRODUCTION
Gluons and quarks are fundamental degrees of freedom of Quantum Chromodynamics (QCD). It is expected that gluons can also serve as building blocks to form hadrons. In the quark model picture, the hadrons made up of valence quarks and valence gluons are usually called hybrids. The hybrid mesons with J P C = 1 −+ are most intriguing since this quantum number is prohibited for qq states of quark model. Up to now there are three experimental candidates for I G J P C = 1 − 1 −+ light hybrid mesons, namely, π 1 (1400) [1], π 1 (1600) [2][3][4] and π 1 (2105) [2] (details can be found in the latest review [5] and the references therein), while lattice QCD studies [6][7][8][9][10][11][12][13] predict that the mass of isovector 1 −+ hybrid meson has a mass around 1.7-2.2 GeV for light quark masses in a range up to the strange quark mass. Very recently, the BESIII collaboration reported the first observation of a I G J P C = 0 + 1 −+ structure η 1 (1855) through the partial wave analysis of the J/ψ → γηη process [14,15]. The resonance parameters of η 1 (1855) are determined to be m η1 = 1855 ± 9 +6 −1 MeV and Γ η1 = 188 ± 18 +3 −8 MeV, and the branching fraction Br(J/ψ → γη 1 (1855) → γηη ) is (2.70 ± 0.41 +0. 16 −0.35 ) × 10 −6 . There have been several phenomenological studies on the properties of η 1 (1855) by assuming it to be an isoscalar light hybrid [16][17][18], a KK 1 (1400) molecular state [19,20], or a tetraquark state [21,22]. As far as the hybrid assignment is concerned, there should be two isoscalar 1 −+ mesons in the flavor SU(3) nonet, and a N f = 2 + 1 lattice QCD study [12] does observe two states of masses around 2.16 GeV and 2.33 GeV, respectively in the 0 + 1 −+ channel (note the light quark mass here corresponds to a pion mass m π ∼ 390 MeV). It is noticed that BESIII also * chenfy@ihep.ac.cn † jiangxiangyu@ihep.ac.cn ‡ cheny@ihep.ac.cn reports a 1 −+ state around 2.2 GeV in the same channel with statistical significance 4.4σ [15]. Since η 1 (1855) is observed in the J/ψ radiative decay, with regard to the possible hybrid assignment, it is desirable to know the production property of the 1 −+ hybrid meson (named as η 1 also) in this process, which will provide important information to the nature of η 1 (1855). This can be investigated in the lattice QCD formalism through the approach similar to the cases of qq mesons [23] and glueballs [24][25][26] in J/ψ radiative decays. The key task is to extract the related electromagnetic multipole form factors from the corresponding three-point functions with a vector current insertion, which involve obviously the annihilation diagrams of the light u, d quarks. Therefore, we adopt the distillation method [27] in the practical calculation, which provides a sophisticated scheme for the operator construction and the computation of all-to-all quark propagators. This paper is organized as follows: Section II presents the details of the numerical calculations of three-point functions, the extraction of the form factors and the interpolation of the form factors to on-shell ones. The discussion of the phenomenological implications of our results can be found in Sec. III. Section IV is a brief summary.

II. NUMERICAL DETAILS
A large statistics is mandatory for the study of J/ψ radiative decay into light hadrons. Our gauge ensemble of N f = 2 degenerate u, d quarks includes 6991 gauge configurations, which are generated on an L 3 ×T = 16 3 × 128 anisotropic lattice with the anisotropy parameter ξ = a s /a t = 5.3 (a s and a t are the spatial and temporal lattice spacing, respectively) [28]. The sea quark mass is tuned to give the pion mass m π ≈ 350 MeV. The parameters of the gauge ensemble are collected in Table I. For the valence charm quark, we adopt the clover fermion action in Ref. [29] and the charm quark mass parameter The partial decay width of J/ψ → γη 1 is governed by the on-shell electromagnetic form factors M 1 (Q 2 = 0) and E 2 (Q 2 = 0) (Q 2 = −p 2 γ ), namely, where α = 1/134 is the fine structure constant at the charm quark mass scale, p γ is the momentum of the final state photon with | p γ | = (m 2 J/ψ − m 2 η1 )/(2m J/ψ ) in the rest frame of J/ψ. These on-shell form factors can be obtained by the Q 2 → 0 interpolation or extrapolation of the form factors M 1 (Q 2 ) and E 2 (Q 2 ), which are defined through the multipole decomposition of the transition matrix elements η 1 (p , λ )|j µ em (0)|J/ψ(p, λ) (see appendix A and also Ref. [30,31]). These matrix elements can be extracted from the following three-point functions where j µ em is the electromagnetic current of quarks, O i η1 ( p, t) and O j J/ψ ( p, t) are the interpolation operators generating η 1 and J/ψ states with a spatial momentum p. Therefore, the major numerical task is to calculate these three-point functions from lattice QCD.
Our lattice setup has the exact SU(2) isospin symmetry. The lattice operator for the isoscalar η 1 takes the form where the chromomagnetic field strength B k is constructed by the proper combination of the gauge covariant spatial derivatives on the lattice [12]. For the operator O i J/ψ we use the conventionalcγ i c-type operator. In order to avoid the complication that the momentum projected operator O i η1 can couple to states with quantum numbers other than 1 −+ [32], the three-point functions in Eq. (2) are calculated practically in the rest frame of η 1 with J/ψ moving at different spatial momenta p. It has been tested that the dispersion relation of J/ψ satisfies FIG. 1. The schematic diagram of the process J/ψ → γη1.
the continuum form very well for all the p modes involved [28]. We only consider the initial state radiation and ignore the case that the photon is emitted from quarks in the final state, so the electromagnetic current j µ em involves charm quarks, namely, j µ em = Z Vc γ µ c (the electric charge of the charm quark Q c = 2 3 e has been absorbed in the prefactor in Eq. (1)). Here Z V is the renormalization constant of the current, since j µ em is not a conserved vector current operator on the lattice. In practice, only the spatial components of j µ em is involved, and its renormalization constant Z s V = 1.118(4) [23] is incorporated implicitly into the expressions in the rest part of this work. Figure 1 illustrates the schematic diagram of Γ iµj after Wick's contraction. It has two separated quark loops, which are actually connected by gluons. The light quark loop on the right hand side can be calculated in the framework of the distillation method. The left part comes from the product of O J/ψ and the current j µ em , namely, (3) which looks very similar to a conventional two-point function of J/ψ and can be calculated independently on each gauge configuration. However, in order for Γ to have good enough signals, the calculation of G µi is highly nontrivial. The conventional momentum source technique turns out to be unfeasible here, because the resulted three-point functions (4) are too noisy even though we have a large gauge ensemble and average over all the time slices τ .
In order to circumvent this difficulty, we calculate G µi in the framework of the distillation method. The distillation method provides a gauge covariant smearing scheme for quark fields, taking the charm quark field where V (t) is the matrix whose columns are eigenvectors of the lattice Laplacian operator −∇ 2 (t) at t (we use N (c) vec = 50 vectors for charm quarks). Therefore, we calculate G µi , whose explicit expression for source time slice at τ = 0 is where S c = cc U is the all-to-all propagator of charm quark for the gauge configuration U and D( p) is a 3L 3 × 3L 3 diagonal matrix with the diagonal matrix elements being δ ij e i p· y ( y labels the column or row indices and i, j = 1, 2, 3 refer to the color indices). Here we apply the

B. Extraction of form factors
When t t 0, the three-point function Γ (3) iµj ( p, 0; t, t ) can be parameterized as where Z η1 ( 0) and Z J/ψ ( p) come from the matrix element which is encoded with the multipole form factors Obviously, in order to extract M iµj ( p), we should know the parameters Z X ( p), m η1 and E J/ψ ( p), which are FIG. 2. Effective mass of π1 (isovector) and η1 (isoscalar), where the fit ranges are [10,30] and [4,14] respectively. The shaded curves illustrate the best-fit values with errors using two-mass-term fits.
actually included in the two-point functions of η 1 and J/ψ, namely, where X stands for J/ψ or η 1 and the source time slice τ is averaged to increase the statistics. The operators O X must be the same as those in the three-point functions Γ iµj , therefore Γ (2) X ( p, t)'s are calculated with the distillation method as well. Since η 1 is set to be at rest, we only calculate Γ (2) η1 ( p, t) at p = 0. The effective mass plot is shown in Fig. 2, where the effective mass of the isovector 1 −+ hybrid state (usually named π 1 ) is also plotted for comparison. The effective mass of η 1 has a much worse signal than that of π 1 due to the inclusion of disconnected diagrams. Through two-mass-term fits in the time range t ∈ [4,14] for η 1 and t ∈ [10,30] for π 1 , the masses are determined to be m π1 = 1.950(28)GeV and m η1 = 2.230(39)GeV, respectively. These results are consistent with those in Ref. [12].
In order for Q 2 to cover the range around Q 2 = 0, the spatial momentum p = 2π Las n of J/ψ is set to run through all possible modes with | n| 2 ≤ 9 based on the m η1 obtained above.
iµj is a smeared operator with N (c) vec = 50, we also generate the perambulators of the valence charm quark with the same N (c) vec to calculate Γ J/ψ ( p) of J/ψ for all the momentum modes involved can be precisely extracted from Γ (2) J/ψ ( p, t) through twomass-term fits. Fig. 3 shows the effective energies E( p, t) (data points) and the fits (colored bands) at different momentum modes n up to | n| 2 = 9.
Along with the calculated two-point functions of J/ψ and η 1 , the matrix element M iµj ( p) is extracted from

E(p, t)/GeV
The effective energies of J/ψ at different spatial modes with | n| 2 ≤ 9. The data points are the numerical results from Γ the ratio function which suppresses the contamination from higher states and should be independent of t and t when ground states dominate. We then make a weighted average value of the function on t to get larger statistics, and take a convention ∆t = t − t .
where ∆ iµj M is the error of the corresponding ratio function, and the weight is to make the average value equal to the least square fit result using a constant. t ∈ [20,40] indicates the "fitting window" in this step. Subsequently, We can extract the form factors M 1 (Q 2 , ∆t) and E 2 (Q 2 , ∆t) from the linear combination of matrix elements M iµj (Q 2 , ∆t) with specific values of i, µ and j. Thus, we can get a similar parameterization for form factors where F i refers to M 1 or E 2 . Note that Q 2 is related to p by iµj ( p, 0; t, t ) is contributed totally by the disconnected quark diagrams, the signal of M iµj ( p, t + ∆t, t ) becomes very noisy when ∆t > ∼ 10 and before a clear plateau appears. Therefore, the resulted M 1 (Q 2 , ∆t) and E 2 (Q 2 , ∆t) have residual time dependence which is absorbed in an additional exponential term in (12). We use this equation as the fitting formula to obtain the value of F i and the corresponding error is acquired from jackknife resampling. Fig. 4 shows the ∆t dependency of M 1 (Q 2 , ∆t) and E 2 (Q 2 , ∆t), whose δm values are listed in Table II.
The fitted parameters, such as M 1 (Q 2 ), E 2 (Q 2 ) and δm are listed in Table II, where one can see that the values of δm at different Q 2 ( n 2 ) are more or less the same value around 1.2-1.3 GeV. This seems a reasonable value. There are quenched lattice QCD calculations of the masses of the first excited 1 −+ strangeonium-like [34] and charmonium-like states [35], which show that the mass differences of the first excited hybrid states and the ground state hybrids are roughly 1.2-1.3 GeV.

C. On-shell form factors and partial decay width
After they are determined at different values of Q 2 , M 1 (Q 2 ) and E 2 (Q 2 ) should be interpolated to the onshell values at Q 2 = 0, which are required to predict the partial decay width using Eq. (1). If a new Lorentz invariant variable is introduced, one can shown that M 1 (Q 2 ) and E 2 (Q 2 ) are proportional to √ Ω, namely, where the form factors G 1 (Q 2 ) and G 2 (Q 2 ) are defined in Eq. (A1) (for details see Appendix A). Obviously, M 1 (Q 2 ) and E 2 (Q 2 ) go to zero when √ Ω → 0. This provides an additional constraint for the Q 2 -interpolation. When putting m = m J/ψ and m = m η1 back to the above expressions, we have Ω(Q 2 ) = m η1 | q| in the rest frame of η 1 . Therefore, it is convenient to introduce a dimensionless function of Q 2 v(Q 2 ) ≡ Ω(Q 2 )/(m J/ψ m η1 ) = | q|/m J/ψ , whose maximum value is v max (Q 2 ) ⇡ 0.49 for the momentumq involved in this study. Note that the form factors G i (Q 2 ) have no singularities when Q 2 > m 2 J/ . They can be expressed as polynomials of Q 2 , and certainly polynomials of v 2 (Q 2 ) (16) where the terms up to O(v 4 (Q 2 )) are kept, since our kinematic configuration that ⌘ 1 is at rest and J/ moves with a momentumq, we have a dimensionless quantity, the velocity of J/ , v(Q 2 ) = p ⌦/(m ⌘1 m J/ ) = |q|/m J/ < 0.494 for the values ofq involved, v 6 max (Q 2 ) ⇠ 1.4% is already much smaller than our statistical errors. Finally, using Eq. (14) we have the interpolation functions for M 1 and E 2 with the constraint F i (Q 2 ) = 0 at v(Q 2 ) = |q|/m J/ = 0 in the rest frame of ⌘ 1 , as suggested by Ref. [36]. Figure 5 shows the Q 2 dependence of M 1 (Q 2 ) and E 2 (Q 2 ) in the region where we are working. The interpolation using Eq. (17) is also illustrated as shaded band in Fig. 5 with the best-fit parameters (the width of the band shows the interpolation error). Thus we get Putting these values into Eq. (1), the partial width is predicted to be (J/ ! ⌘ 1 ) = 2.04(61) eV (19) (using the ⌘ 1 mass m ⌘1 = 2.230(39) GeV). Note that the form factors in Eq. (18) are obtained by assuming ⌘ 1 to be a stable particle, while it must be a resonance in principle. For a resonance R of parameters (m R , R ), a more systematic approach to derive the form factor f R (Q 2 ) from lattice QCD has been proposed in Refs. [37][38][39][40][41][42] where the finite volume correction are thoroughly discussed for the 1+J ! 2 type transitions with J being a local current, especially for the case that a resonance can appear in the final two hadron system. However, this approach is unfeasible yet for the processes J/ ! + lighthadron(s) that take place solely through quark annihilation diagrams, because the low precision of (3) iµj cannot a↵ord that sophisticated treatment. Fortunately, some examples [37,41,42] indicate that the finite volume correction to the form factors of a narrow resonance R is O( R /m R ) when R is treated as a stable particle. If ⌘1 /m ⌘1 in the N f = 2 case is similar to or even smaller than that of ⌘ 1 (1855), the form factors in Eq. (18) may be taken as approximations for those of the resonant ⌘ 1 with regard to their large statistical uncertainties of roughly 15%-20%. whose maximum value is v max (Q 2 ) ≈ 0.49 for the momentum q involved in this study. Note that the form factors G i (Q 2 ) have no singularities when Q 2 > −m 2 J/ψ . They can be expressed as polynomials of Q 2 , and certainly polynomials of v 2 (Q 2 ) (16) where the terms up to O(v 4 (Q 2 )) are kept, since our kinematic configuration that η 1 is at rest and J/ψ moves with a momentum q, we have a dimensionless quantity, the velocity of J/ψ, v(Q 2 ) = √ Ω/(m η1 m J/ψ ) = | q|/m J/ψ < 0.494 for the values of q involved, v 6 max (Q 2 ) ∼ 1.4% is already much smaller than our statistical errors. Finally, using Eq. (14) we have the interpolation functions for M 1 and E 2 with the constraint F i (Q 2 ) = 0 at v(Q 2 ) = | q|/m J/ψ = 0 in the rest frame of η 1 , as suggested by Ref. [36]. Figure 5 shows the Q 2 dependence of M 1 (Q 2 ) and E 2 (Q 2 ) in the region where we are working. The interpolation using Eq. (17) is also illustrated as shaded band in Fig. 5 with the best-fit parameters (the width of the band shows the interpolation error). Thus we get Putting these values into Eq. (1), the partial width is predicted to be Γ(J/ψ → γη 1 ) = 2.04(61) eV (19) (using the η 1 mass m η1 = 2.230(39) GeV). Note that the form factors in Eq. (18) are obtained by assuming η 1 to be a stable particle, while it must be a resonance in principle. For a resonance R of parameters (m R , Γ R ), a more systematic approach to derive the form factor f R (Q 2 ) from lattice QCD has been proposed in Refs. [37][38][39][40][41][42] where the finite volume correction are thoroughly discussed for the 1+J → 2 type transitions with J being a local current, especially for the case that a resonance can appear in the final two hadron system. However, this approach is unfeasible yet for the processes J/ψ → γ + lighthadron(s) that take place solely through quark annihilation diagrams, because the low precision of Γ  [37,41,42] indicate that the finite volume correction to the form factors of a narrow resonance R is O(Γ R /m R ) when R is treated as a stable particle. If Γ η1 /m η1 in the N f = 2 case is similar to or even smaller than that of η 1 (1855), the form factors in Eq. (18) may be taken as approximations for those of the resonant η 1 with regard to their large statistical uncertainties of roughly 15%-20%.

III. DISCUSSION
Although obtained for N f = 2, the form factors in Eq. (18) can be applied to the discussion of the physical SU(3) case. In J/ψ radiative decays, the final state light hadron (η 1 here) is produced by the gluons from the cc annihilation, and thereby must be a flavor singlet (isoscalar for N f = 2 and SU (3) singlet for N f = 3). If the flavor wave function of the light hadron is properly normalized, the underlying gluonic dynamics is usually independent of N f except for the U A (1) anomaly relevant interaction. In this sense, the form factors in Eq. (18) On the other hand, the masses of η (l,h) 1 can be different from m η1 in this study, we should consider the correction factor due to the mass mismatch. According to Eq. (14), one has Since the form factors G i (Q 2 ) are functions of Q 2 and are regular around Q 2 = 0, it is expected the form factors G i (0) for i = 1, 2 are insensitive to m η1 in the range m η1 ∼ 2 GeV. For the case of this study, G 2 1 (0) is a few times larger than G 2 2 (0), such that from Eq. (1), the m η1 dependence is approximately Γ ∝ | q| 3 /m 2 η1 . Thus one has the following partial widths, | pγ (η1)| 3 is the compensating kinematic factor due to the mass mismatch of η 1 and η (x) 1 . As for the ηη decay mode where η 1 (1855) is observed, since it must be a flavor octet, the flavor SU(3) symmetry implies the decay η = sin θ ηη |H I |η (8) where g is the effective coupling, is the polarization vector of η (x) 1 and k (x) is the momentum of ηη in the η (x) decay. Thus we obtain the ratio (24) which is free from θ but depends solely on the masses and widths of η is not too large, the kinematic factor in the above The lattice QCD study in Ref. [12] observes η or equivalently the admixtures of ss and nn = (uū + dd)/ √ 2 through a mixing angle α, If the flavor wave functions of η one can easily show that θ is related to α by θ = α−54.7 • . This convention for the mixing angle α is the same as that in Ref. [12] where α is determined to be roughly α = 22.7(2.1) • (averaged over the values on the three lattices involved), such that one has θ ≈ −32.0(2.1) • . This indicates a large mixing of η couples strongly to ηη , namely, the effective coupling in Eq. (23) is roughly g = 5.0(1.0) (note the effective coupling g ρππ ≈ 6.0 for the decay process ρ → ππ). Although the possible enhancement by QCD U A (1) anomaly [16,17], this is really a large coupling and should be understood when comparing with the significantly small coupling of its isovector partner π 1 to η π, which is expected by phenomenological studies [44,45] and estimated by lattice QCD calculations [11,13].

IV. SUMMARY
Based on a large gauge ensemble of N f = 2 dynamical quarks at m π ≈ 350 MeV, we perform the first theoretical calculation of Γ(J/ψ → γη 1 ) where η 1 is the light flavor singlet 1 −+ hybrid. The related three-point functions are contributed totally from disconnected quark diagrams, which are dealt with using the distillation method. The on-shell electromagnetic form factors are determined to be M 1 (0) = −4.73(74) MeV and E 2 (0) = 1.18 MeV, which give Γ(J/ψ → γη 1 ) = 2.04(61) eV for m η1 = 2.23(4) GeV. These results are applicable to discuss the production rates of the two mass eigenstates η in the SU(3) case, if the singlet-octet mixing angle is known. As for η 1 (1855) observed by BESIII, its hybrid assignment depends strongly on the existence of its mass partner. It should be emphasized that the ratio of the branching fractions Br(J/ψ → γη (l,h) 1 → γηη ) is inversely proportional to the ratio of the total widths of η (l,h) 1 . This can be used as one of the criteria to identify η (l,h) experimentally. If η 1 (1855) is a hybrid for sure, our results and the mixing angle θ determined in Ref. [12] indicate that the coupling of the octet 1 −+ hybrid η (8) 1 to ηη is very large. This is interesting and worthy of an investigation in depth. Throughout our calculation, η 1 is tentatively viewed as a stable particle. This surely introduce theoretical uncertainties which cannot be accessed in the present stage, but should be explored in future works. Nevertheless, this study provides the first valuable theoretical predictions for this intriguing topic from lattice QCD.  [46] and QUDA library [47,48] are acknowledged. The computations were performed on the HPC clusters at Institute of High Energy Physics (Beijing) and China Spallation Neutron Source (Dongguan), and the ORISE computing environment.