Three Photon Decay of $J/\psi$ from Lattice QCD

The three photon decay rate of $J/\psi$ is studied using a $N_f=2+1+1$ twisted mass gauge ensemble with $60$ configurations and the results are compared with the existing experiments by CLEOc and BESIII. Only the correlation function that is directly related to the physical decay width is computed, with all polarization of the initial and final states summed over. We obtain the result for such rare decay $\mathcal{B}(J/\psi\rightarrow 3\gamma)=(3.62 \pm 0.04\pm 1.50)\times 10^{-5}$ where the first error is statistical and the second is the estimated systematic error. This is in the right ballpark when compared with both CLEOc and BESIII. Furthermore, based on our lattice result, an interesting enhancement structure is discovered in the Dalitz plot which can be verified in future experiments with more statistics.

The modern tools of treating the quarkonium spectrum, decay and production are nonrelativistic QCD (NRQCD) and factorization [6], where the decay rate of J/ψ → 3γ are expanded by the lowest order NRQCD J/ψ-to-vacuum matrix element plus relativistic corrections. However, it has been known that perturbative calculations for this process suffers from both large and negative radiative and relativistic corrections [7]. Both the inconsistency between theory and experiments and the divergence puzzle appearing in higher order radiative corrections [8] indicate that the NRQCD might break down for predicting the decay rate of J/ψ → 3γ. Although three decades have passed, no much better understand-ing of the process J/ψ → 3γ developed by NRQCD than the potential model in early 1980s [9][10][11].
Therefore, it is of great significance to seek for new methods. In this letter, we propose to use lattice QCD as such an alternative and present the first exploratory computation of J/ψ → 3γ decay width using an ensemble of twisted mass gauge field configurations. In lattice QCD calculations, one normally evaluates the matrix element of an interested interpolating operator with correct quantum numbers between two hadronic states. Though photon is not an eigenstate of QCD, regarding the photon as a superposition of QCD eigenstate and adopting the corresponding eletromagnetic current J µ em as the effective photon interpolating operators has been put forward [12]. This amounts to treat QED perturbatively while keeping the QCD part nonperturbative. This method has been widely used in photon structure functions [13], radiative transition [14] and two-photon decays in charmonia [15,16].
Method -We start by expressing the amplitude of J/ψ → 3γ in terms of the appropriate four-point function using Lehmann-Symanzik-Zimmermann reduction formula, integrating out the photon fields perturbatively and continuing the resulting expression to Euclidean space analytically. This process introduces the photon virtualities Q 2 i = |q i | 2 − ω 2 i (for more details, see Ref. [15]), we then arrive at the following final result for the four-point function that is relevant for the process J/ψ → 3γ, Here the four polarization vectors: µ , ν , ρ and α correspond to the three final photons and the initial J/ψ particle, respectively, with the polarizations labelled by λ 1 , λ 2 , λ 3 and λ 0 . The analytic continuation from Minkowski to Euclidean space here works out as long as the virtualities of three photons are not too timelike to produce on-shell vector hadrons. More specifically, where M V is mass of the lightest vector meson. For simplicity, we have used the local current j µ (x) =c(x)γ µ c(x) for the charm quark which can be renormalized by a multiplicative factor Z V , the value of which is available in the literature, namely, Z V = 0.6095(03) [17]. The correlation functions appearing in the above equation can then be evaluated in lattice QCD in terms of quark propagators. In this exploratory calculation, we have neglected the disconnected diagrams therefore only charm propagators are needed. We now denote the matrix ele-  µνρσ |M µνρσ | 2 which will be called Tfunction in the following. As we will see, T -function actually represents a distribution of physical partial decay width in terms of a pair of kinematic variables. Each M µνρα can be computed on the lattice using the fact that M is independent of the time t, as long as |t f − t| is large enough.
In real simulations, plateau behaviours are searched for to extract the values of M in various cases. The current that couples to the first photon is fixed at t i , other two are placed at t and t, respectively and J/ψ meson is placed at t f (as shows in Figure. 1). The integrals in Eq. (1) are also replaced by a trapezoidal summation.
In conventional lattice computations, for example in the decay of η c → γγ etc., the hadronic matrix element such as M µνρα is further decomposed into various form factors (which are functions of the virtualities Q 2 i ) according to its Lorentz structure. By fitting the matrix element at different Q 2 i with a particular functional form, one arrives at the complete off-shell form factors and finally, by setting all virtualities to the onshell values, namely Q 2 i = 0, one obtains the final physical decay rate. In our case, the form factor decomposition is way too complicated. In the study of 3γ decays of Z and the positronium, people have worked out such decompositions in perturbation theory [18] [19]. However, there is no guarantee that these perturbative decompositions will also work in QCD. Therefore, we will proceed in another way. We will be satisfied with the physical decay width only. That is to say, we will be only interested in the on-shell matrix element. Thus, we can perform the summation over polarizations of the initial and final particles first, and only the on-shell matrix element will be computed on the lattice. Due to Ward-identities of the currents, the summation over polarizations of the photons yields the Minkowski metric, e.g. λ1 ( µ (q 1 , λ 1 ) * µ (q 1 , λ 1 ) ⇒ η µµ . The summation over the initial polarization of J/ψ yields the same if we take the rest frame of the particle. Therefore, we have λi |M| 2 = µνρα |M µνρα | 2 . In our actual simulations, we sum over all physical polarizations (altogether 24 possibilities) of M.
The decay width of J/ψ → 3γ in J/ψ center of mass frame can be expressed as, where x, y are two dimensionless variables in the range It is easily checked that they fall into the right-upper triangle of the unit square in the xy-plane, i.e. satisfying : x ∈ [0, 1], y ∈ [1 − x, 1]. In the continuum, the on-shell three-body decay pattern are normally parameterized by the so-called Dalitz plots, which can be obtained from the T -function T (x, y).
However, due to the discreteness of the momenta on the finite lattice, it is impossible to exactly impose the on-shell condition for all particles, making the on-shell quantity T (x, y) not directly accessible. Instead, the on-shell conditions for the particles can be done as follows: We first put the J/ψ particle and at least one final photon on shell, keeping the other two photons as close to on-shell as possible by adjusting their three momenta. It is found that this still introduce some non-vanishing virtualities to the other photons. With these non-vanishing but small virtualities, the matrix element can be computed directly on the lattice, the norm of which we denote as T (x, y, Q 2 1 , Q 2 2 , Q 2 3 ). This differs from the T -function only because of the fact that some of the photons are still not on-shell. We then try to estimate the on-shell quantity, the T -function T (x, y), by the following fitting formula, for |Q 2 i | 1 where everything is measured in lattice units. We expect such behavior since the final three photons are identical particles.  (2) and (3) constitute the central part of this letter. As pointed out already, different from the conventional method, we have intentionally avoided the amplitude parameterization for J/ψ → 3γ, though it has the similar structure as Z → 3γ [18] and the positronium to 3γ decay [19]. Because all the form factors introduced in the amplitude parameterizations are scalar functions of three-photon momenta, permutations of these momenta then lead to more form factors, rendering the computation of all of these form factors too costly. In the case of three-body decay, what is really measured in the experiments are the so-called Dalitz plots of the final states. Dalitz plot represents the distribution of the partial decay width in two independent kinematic variables. In the case of three-photon decay of J/ψ, this is taken to be the largest and the smallest two-photon invariant mass values, denoted as M (γγ) lg and M (γγ) sm , respectively, among three combinations for the final photons. We will call them the Dalitz variables in the following. These two Dalitz variables are directly related to the kinematic variables (x, y) that we introduced. To be more specific, we have, where the upper/lower line on the right corresponds to the case of M (γγ) lg /M (γγ) sm , respectively. Thus, the Dalitz plot for the three-body decay is directly related to the on-shell T -function T (x, y) that we aim to compute on the lattice.
Simulations And Results -Our lattice calculation is performed using a N f = 2 + 1 + 1 flavor twisted mass gauge filed ensemble with 60 configurations at coupling β = 1.95. The size of the lattice is 32 3 × 64 with lattice spacing a ≈ 0.0823 fm and light quark mass that corresponds to physical pion mass of m π ∼ 372 MeV which was generated by Extended Twisted Mass Collaboration (ETMC) [17]. The conventional sequential method has been adopted to calculate the four-point function. Two sequential sources are placed close to J/ψ meson, and the contraction is performed on the furthest current. After the integration (summation) of time slice t i and t , the matrix element M µνρα , being a function of time slice t, can be obtained on the lattice. We have chosen 5 sets of suitable photon virtualities Q 2 i and momenta which results in 24 off-shell four-point functions M µνρα for each set. In Fig. 2, typical plateau behaviors for the four-point func-tion M µνρα are shown in the case of µνρα = 1111. The data points with errors are the results from the simulation and the errors are estimated using jackknife method. Other cases are similar. It is seen that with 60 configurations, the signal is reasonable and the statistical errors are under control.
As is seen from Fig. 2, although only 5 sets of photon three-momenta are considered, the physical amplitude J/ψ(p, λ 0 )|γ(q 1 , λ 1 )γ(q 2 , λ 2 )γ(q 3 , λ 3 ) is invariant under the photon exchange (q i , λ i ) ↔ (q j , λ j ), so we finally obtain the off-shell T -function, i.e. T (x, y, Q 2 1 , Q 2 2 , Q 2 3 ), at a total of 27 points of (x, y) in the xy-plane. For each of these, the on-shell function T (x, y) can be extracted by performing a correlated fit using Eq. (3). This is shown in FIG. 3. The correlation among different Q 2 i points are estimated using bootstrap.
For three-body decays, the decay width is a quantity that both experimentalists and theorists are interested in. Most of time, however, the total decay width itself is not directly measurable for experimentalists. They first obtain the so-called Dalitz plot, which is a distribution of the decay width in the plane of two kinematic variables. In the case of J/ψ → 3γ, these are exactly the Dalitz variables M (γγ) lg and M (γγ) sm that we mentioned, which are related to the (x, y) kinematic variables via Eq. (4). Being first-hand data obtained in experiments, Dalitz plot plays a key role for a three-body finial state. As is well-known, bands that appear in the Dalitz plot indicate that there is an intermediate two-body state. Thus, nonuniformity in the Dalitz plot can offer immediate information on the cross section |M| 2 . As we will illustrate below, this can be related to the on-shell T -function that we compute on the lattice.
It is seen that, the relation between the Dalitz variables and the pair (x, y) as indicated in Eq. (4), maps the upper right triangular region of unit square in the (x, y) plane onto a corresponding region in (M (γγ) sm , M (γγ) lg ) plane in the Dalitz plot. The shape of the region in the Dalitz variables is not regular, with a rather sharp tip extending towards the right lower corner in the (M (γγ) sm , M (γγ) lg ) plane. This is exactly what is measured in the experiments, see e.g. Fig.1 (d) in Ref. [5]. As we have already obtained T (x, y) at 27 points in the (x, y) plane, we utilize a cubic spline function to interpolate these points. The surface of this resulting interpolating function T (int) (x, y) is illustrated in the top half of Fig. 4, together with the original data points shown in red. In the lower half of the same figure, we illustrate the mapping from (x, y) plane to the Dalitz variables plot as suggested in Eq. (4). On the left is the contour plot of the interpolated function T (int) (x, y) while on the right is the corresponding one in Dalitz variables. To further obtain the total decay width for the process, one needs to either integrate the function T (x, y) in the (x, y) plane, or doing the corresponding integration in the Dalitz variables.
Before we integrate the function to give the final decay rate, let us make the following comments: i) We have computed the connect diagram as shown in Fig. 1. This diagram can in principle also include the physical process J/ψ → γη c → γγγ as well. Therefore, in order to make comparison with the experiments, we need to remove such contributions from our lattice data. It is easily verified that this corresponds to the three corners of the triangle in the (x, y) plane. In the experiments, these are also the region where the major background comes in. To remove these contributions, we need to make definite cuts as the experimentalists did, see e.g. Ref. [4,5]. For example, if we cut the three corners by the condition M (γγ) lg < 2.9GeV, this results in a deduction of 0.12eV in the final result for Γ(J/ψ → 3γ) quoted in Eq. (5) below.
iii) An equal-energy enhancement structure of three photons for J/ψ → 3γ is observed. There is an apparent enhancement when the energies of three photons are approximately equal for the T -function in the (x, y) plane, located close to the center of the triangle. Transform into the Dalitz variables, such an enhancement appears close to the tip. Although the enhancement is not that dramatic, we expect future experiments would observe this if enough statistics of J/ψ are collected.
Now we proceed to integrate T -function surface, i.e. the function T (int) (x, y), over the physical region and finally arrive at the value for the decay width of J/ψ → 3γ, with the contribution of 0.12eV from γη c already subtracted . Here the first error 0.04 only accounts for statistical errors, part of which is estimated by resampling the original T (x i , y i ) data using bootstrap method. Also included here is the on-shell fitting process as suggested in Eq. (3). The above mentioned errors are added up by quadrature yielding 0.04. The second error 1.40 represents our rough estimate of the systematic errors. There are two sources of the systematic errors that we estimated: there is a part with magnitude 0.06 arising from the cubic spline interpolation process; the rest arises from the twisted-mass term, which comes from the lattice discretization. This is estimated by swapping the three-momentum q 1 and q 3 of the two final photons. Physically, this should result in the same matrix element M. However, since twisted mass fermion endows lattice artifacts at the order of O(a 2 ), this results in a difference in the matrix element hence a different T -function. Numerically, this difference is found to be as large as 38% for the T -function, which is somewhat surprising. Needless to say, there are more systematic errors that need to be taken into account. Hopefully, these can be considered in future more systematic lattice studies.
The branching fraction, if the uncertainty of J/ψ total width ignored, is given by B(J/ψ → 3γ) = (3.62 ± 0.04±1.50)×10 −5 , to be compared with the BESIII result (1.13 ± 0.18 ± 0.2) × 10 −5 [5]. Our lattice result is larger than both the results of CLEOc and BESIII Collaboration, however, bearing in mind that our systematic error is only an estimation, it is still reassuring to see that our preliminary lattice result is in the right ballpark. Needless to say, in the future, more systematic lattice studies are much welcome.
Discussions -As we have said, the central value of our lattice result for the decay rate Γ(J/ψ → 3γ) turns out to be larger than the experimental results. More over, we have also discovered a enhancement structure in the corresponding Dalitz plot. As we argue below, we suspect that such a structure has not been observed in the experimental Dalitz plot so far is likely due to the lack of statistics. In what follows, we briefly discuss the cause of these discrepancies. First, concerning the total decay rate, the events of J/ψ decay collected in experiments include not only the direct decay mode J/ψ → 3γ, but also other modes. For example, the prominent contributions J/ψ → γ/π 0 /η/η /η c → 3γ. To avoid these effects, events lie in corresponding energy bands are rejected in the experiment. Inevitably, some events from the direct decay in such bands are also rejected, and these contributions may not be small, as is seen from the corresponding band covering different areas in Figure. 4 for the Dalitz plot. By contrast, the process in Figure. 1 is cleaner.
Second, let us move on to the comparison of the Dalitz plot, which is a more detailed comparison of the lattice result with the experiment. For such a purpose, we define a normalized T -function distribution densityT (x, y) as which can be viewed as probability density in the (x, y) plane. The corresponding Dalitz plot can also be generated by drawing random samples using this probability distribution. In Fig. 5 Therefore, we believe that it is the lack of statistics that leads to the disappearance of equal-energy enhancement region in the experiments. If such a structure is for real, we expect BESIII would be able to observe it in the future with 1.31 × 10 9 J/ψ events already collected. Third, it is observed that the statistical errors in our lattice calculation are negligible when compared with the systematic ones. Therefore, in the future, it becomes crucial to reduce the systematic errors so as to have a closer comparison. Since we have witnessed a surprisingly large error that arises due the violation of parity of twisted-mass fermions, it is natural to implement a similar computation using other type of fermions.
Finally, the lattice calculation presented in this letter can promote other similar lattice computations. By summing over final and initial state polarizations, one can compute the distribution of the partial decay width in the corresponding Dalitz plot, which can be compared directly with the experiments. In principle, we could also keep the information of the initial polarization of the J/ψ particle, which would give us the decay rate for polarized J/ψ samples. Of course, one could also contemplate to generalize the current framework to other hadronic decays with three particles in the final state.
Conslusions -Using lattice QCD, an exploratory calculation of rare decay rate Γ(J/ψ → 3γ) is presented. We obtain the branching fraction B(J/ψ → 3γ) = (3.62 ± 0.04 ± 1.50) × 10 −5 where the first error is statistical and the second being an estimate of the systematics. The central value turns out to be larger than the two existing experimental ones from CLEOc and BESIII. It is still reassuring to see that our lattice result is in the right ballpark. We also briefly discuss the possible causes of the discrepancy. It is also seen that the current lattice calculation is dominated by systematic effects, with statistical errors rather well-controlled.
Instead of parameterizing the matrix relevant matrix element with form factors, we evaluate the squared matrix element which is directly related to the physical decay width. We could even obtain the distribution of the decay width using two kinematic variables that is related to the Dalitz plot in the experiments. For the decay J/ψ → 3γ we have also found a equal-energy enhancement region in the Dalitz plot which could be tested by BESIII Collaboration in the future, with their largest J/ψ sample in the world.