Some hadronic parameters of charmonia in $\boldsymbol{N_{\text{f}}=2}$ lattice QCD

The phenomenology of leptonic decays of quarkonia holds many interesting features: for instance, it can establish constraints on scenarios beyond the Standard Model with the Higgs sector enriched by a light CP-odd state. In the following paper, we report on a two-flavor lattice QCD study of the $\eta_c$ and $J/\psi$ decay constants, $f_{\eta_c}$ and $f_{J\psi}$. We also examine some properties of the first radial excitation $\eta_c(2S)$ and $\psi(2S)$.


I. INTRODUCTION
The discovery at LHC of the Higgs boson with a mass of 125.09(24) GeV [1] has been a major milestone in the history of Standard Model (SM) tests: the spontaneous breaking of electroweak symmetry generates masses of charged leptons, quarks and weak bosons. A well-known issue with the SM Higgs is that the quartic term in the Higgs Lagrangian induces for the Higgs mass m H a quadratic divergence with the hard scale of the theory: it is related to the so-called hierarchy problem. Several scenarios beyond the SM are proposed to fix that theoretical caveat. Minimal extensions of the Higgs sector contain two complex scalar isodoublets Φ 1,2 which, after the spontaneous breaking of the electroweak symmetry, lead to 2 charged particles H ± , 2 CP-even particles h (SM-like Higgs) and H, and 1 CP-odd particle A. In that class of scenarios, quarks are coupled to the CP-odd Higgs through a pseudoscalar current. Those extensions of the Higgs sector have interesting phenomenological implications, especially as far as pseudoscalar quarkonia are concerned. For example, their leptonic decay is highly suppressed in the SM because it occurs via quantum loops but it can be reinforced by the new tree-level contribution involving the CP-odd Higgs boson, in particular in the region of parameter space where the new boson is light (10 GeV m A 100 GeV) and where the ratio of vacuum expectation values tan β is small (tan β < 10) [2,3]. Any enhanced observation with respect to the SM expectation would be indeed a clear signal of New Physics. Let us finally note that the hadronic inputs, which constrain the CP-odd Higgs coupling to heavy quarks through processes involving quarkonia, are the decay constants f ηc and f η b .
This paper reports an estimate of hadronic parameters in the charmonia sector using lattice QCD with N f = 2 dynamical quarks: namely, the pseudoscalar decay constant f ηc , because of its phenomenological importance, but also the vector decay constant f J/ψ as well as the ratio of masses m ηc(2S) /m ηc and m ψ(2S) /m J/ψ . The two latter quantities are very well measured by experiments and their estimation has helped us to understand how much our analysis method can address the systematic effects, on quarkonia physics, coming from the lattice ensembles we have considered. We present also our findings for the following ratios of decay constants: f ηc(2S) /f ηc and f ψ(2S) /f J/ψ . This work is an intermediary step before going to the bottom sector, the final target of our program, because it is more promising for the phenomenology of extended Higgs sectors. An extensive study of the spectoscopy of charmonia has been done in [4,5] while only two lattice estimates of f ηc and f J/ψ are available so far at N f = 2 [6] and N f = 2 + 1 [7,8] This study has been performed using a subset of the CLS ensembles. These ensembles were generated with N f = 2 nonperturbatively O(a)-improved Wilson-Clover fermions [9,10] and the plaquette gauge action [11] for gluon fields, by using either the DD-HMC algorithm [12][13][14][15] or the MP-HMC algorithm [16]. We collect in Table I our simulation parameters. Two lattice spacings a β=5.5 = 0.04831 (38) fm and a β=5.3 = 0.06531(60) fm, resulting from a fit in the chiral sector [17], are considered. We have taken simulations with pion masses in the range [190 , 440] MeV. The charm quark mass has been tuned after a linear interpolation of m 2 Ds in 1/κ c at its physical value [18], after the fixing of the strange quark mass [19]. The statistical error on raw data is estimated from the jackknife procedure: two successive measurements are sufficiently separated in trajectories along the Monte-Carlo history to neglect autocorrelation effects. Moreover, statistical errors on quantities extrapolated at the physical point are computed as follows. Inspired by the bootstrap prescription, we perform a large set of N event fits of vectors of data whose dimension is the number of CLS ensembles used in our analysis (i.e. n = 6) and where each component i of those vectors is filled with an element randomly chosen among the N bins (i) binned data per ensemble. The variance over the distribution of those N event fit results, obtained with such "random" inputs, is then an estimator of the final statistical error. Finally, we have computed quark propagators through two-point correlation functions using stochastic sources which are different from zero in a single timeslice that changes randomly for each measurement. We have also applied spin dilution and the one-end trick to reduce the stochastic noise [20,21].

GEVP discussion
The two-point correlation functions under investigation read where V is the spatial volume of the lattice, ... is the expectation value over gauge configurations, and the interpolating fieldscΓc can be non-local. As a preparatory work, different possibilities were explored to find the best basis of operators, combining levels of Gaussian smearing, interpolating fields with a covariant derivativecΓ( γ · ∇)c and operators that are odd under time parity. Solving the Generalized Eigenvalue Problem (GEVP) [22,23] is a key point in this analysis. Looking at the literature we have noticed that people tried to mix together the operatorscΓc andcγ 0 Γc in a unique GEVP system [4,6]. In our point of view, this approach raises some questions: let us take the example of the interpolating fields {P =cγ 5 c ; A 0 =cγ 0 γ 5 c}. The asymptotic behaviours of the 2-pt correlation functions defined with these interpolating fields read The matrix of 2 × 2 correlators of the GEVP is then In the general case, the spectral decomposition of C ij (t) is The dual vector u n to Z s is defined by However, in the case of mixing T-odd and T-even operators in a single GEVP, the D's do depend on i and j: the previous formula for λ n (t, t 0 ) is then not correct. Hence, approximating every correlator by sums of exponentials forward in time, together with the assumption that the D ij are independent of i and j, may face caveats. A toy model with 3 states in the spectrum helps to understand this issue: In our numerical application, we have chosen T = 64, t 0 = 3 and compared 2 × 2 and 3 × 3 subsystems: results can be seen in Fig.1. Our observation is that until t = T /4, neglecting the time-backward contribution in the correlation function has no effect. Based on this study, we can affirm that the method is safe for the ground state and the first excitation. In contrast, one may wonder what would happen with a dense spectrum when the energy of the third or a higher excited state is extracted: in that case there will be a competition between the contamination from higher states, which are not properly isolated by the finite GEVP system, and the omission of the backward in time contribution to the generalised eigenvalue under study, even at times t < 1 fm.

Interpolating field basis
Building a basis of operators with can have advantages and this was already explored [24,25]. But there are sometimes bad surprises. . . a good example being the pseudoscalar-pseudoscalar correlator defined by whose behaviour will be anticipated using the quark model method.
In the formalism of quark models, the c quark field reads being developped into a large and a small component. Notice that the c quark field is a spinor of type u while thec antiquark field is of type v, so at this stage we need to expressv in terms of u. Defining the charge conjugation operator by C = −iγ 0 γ 2 and using the Dirac representation for the γ matrices, we can write Thec antiquark field is thenc Assuming the charmonium at rest, we have Then ordinary Pauli matrices algebra leads tō The previous approach is very naive with respect to quantum field theory, in particular the approximation D i c 1 ∼ p i c 1 , but as a conclusion one can see that interpolating fieldsc γ 5 ( γ· ∇)c potentially give very noisy correlators. And, indeed, we have found a numerical cancellation between the "diagonal" contribution and the "off-diagonal" contribution resulting in a very noisy correlator C(t) which is compatible with zero.

Summary
To summarise, we have considered four Gaussian smearing levels for the quark field c, including no smearing, to build 4 × 4 matrix of correlators without any covariant derivative nor operator of the π 2 or ρ 2 kind [24], from which we have also extracted the O(a) improved hadronic quantities we have examined. Solving the GEVP for the pseudoscalar-pseudoscalar and vector-vector matrices of correlators we obtain the correlators that will have the largest overlap with the n th excited state as follows: and their symmetric counterpart with the exchange of operators at the source and at the sink. The quark bilinears are P =cγ 5 c, A 0 =cγ 0 γ 5 c, V k =cγ k c and T k0 =cγ k γ 0 c. Moreover, in those expressions, the label L refers to a local interpolating field while sums over i and j run over the four Gaussian smearing levels. The projected correlators have the following asymptotic behaviour

Considering the O(a) improved operators
where the lattice derivatives are defined by the f ηc decay constants is extracted in the following way: while the f J/ψ decay constant is obtained with: The R superscript denotes the renormalized improved operators. The renormalization constants Z A and Z V have been non perturbatively measured in [26,27]. We have also used non-perturbative results and perturbative formulae from [28][29][30] for the improvement coefficients c A , c V , b A and b V , and the matching coefficient Z between the quark mass m c defined through the axial Ward Identity and its counterpart defined through the vector Ward Identity The decay constants of the considered radial excited states are given by

B. Analysis
Let us first consider the mass and the decay constant of the ground states. Since the fluctuations with time are large, we have decided to use generalized eigenvectors at fixed time t fix , v P (V ) 1 (t fix , t 0 ), in order to perform the corresponding projection. In practice we have chosen t fix /a = t 0 /a + 1 but we have checked that the results do not depend on t fix . However, for the excited states, the formulae written above apply. Although the fluctuations are even larger than for the ground states, the correlatorsC 2 A0P andC 2 V V are qualitatively well fitted by a single contribution, which is not true at all if the projection on v P (V ) 2 (t fix , t 0 ) is used: effective masses extracted with the latter projection show a very big slope with time. For the ground states, the time range [t min , t max ] used to fit the projected correlators is set so that the statistical error on the effective mass δm stat (t min ) is larger enough than the systematic error ∆m sys (t min ) ≡ exp[−∆t min ] with ∆ = E 4 − E 1 ∼ 2 GeV. That guesstimate is based on our 4 × 4 GEVP analysis, though we do not want to claim here that we really control the energy-level of the third excited state. Actually our criterion is rather δm stat (t min ) > 4∆m sys (t min ) to be more conservative. On the other side, t max is set after a qualitative inspection of the data which count for the plateau determination. Since ∆ (P ) ∼ ∆ (V ) , fit intervals are identical for pseudoscalar and vector charmonia. For the first excited states, the time range has been set by looking at effective masses and where the plateaus start and end, including statistical uncertainty. Finally, t 0 = 3a at β = 5.3 and t 0 = 5a at β = 5.5: the landscape is unchanged by increasing t 0 . We have collected in Table II of the Appendix the raw data we have obtained in our analysis. The extrapolation to the physical point of the quantities we have measured uses a linear expression in m 2 π with inserted cut-off effects in a 2 : X(a, m π ) = X 0 + X 1 m 2 π + X 2 (a/a β=5.3 ) 2 We remind that κ c have been tuned at every κ sea so that m Ds (κ s , κ c , κ sea ) = m phys Ds (it was necessary to tune the strange quark mass κ s , though it is not a relevant quantity for the study discussed here). The Fig.2 shows the effective masses of the η c , η c (2S), J/ψ and ψ(2S) mesons for the set F7 which are obtained by the following expressions The drawn lines correspond to our plateaus for the fitted masses. One can observe that our plateaus are large for the ground states but they are unfortunately of a much worse quality for the radial excitations. The latter data also show large fluctuations with time, for which we did not find any satisfying explanation. If t fix were used in the correlators projection procedure, the price to pay would be a strong contamination by other states inC 2 A0P andC 2 V V other than the ones of interest. We show in Fig.3  together with the experimental determination of the J/ψ mass and width, and setting α em (m 2 c ) = 1/134 [31], one gets f "exp" J/ψ = 407(6) MeV We collect the various lattice QCD results and the phenomenological estimate of f J/ψ in Fig.5. Unfortunately the situation is less favorable with the excited states. Fig.6 shows the extrapolation to the continuum limit of the ratios m ηc(2S) /m ηc and m ψ(2S) /m J/ψ , compared to the experimental values. Since the cut-off effects are small, of the order of 5%, there is no hope to have points in the continuum limit significantly lower than our lattice data. We get m ηc(2S) /m ηc = 1.281 (7) (m ηc(2S) /m ηc ) exp = 1.220 In [32], lattice results were much closer to the experiment but the slope in a 2 was probably overrestimated from the coarsest lattice point: points at lattice spacings similar to those used in our work are compatible with our data. The situation is even more confusing for the ratios of decay constants, as shown in Fig.7 f ηc(2S) /f ηc = 0.81(8) < 1 while We have not found any explanation for such a large spin breaking effect. Projected correlatorsC 2 in the pseudoscalar and the vector sector show the same quality fit with similar fluctuations. We also have performed a global fit of individual correlators C V (i) V (j) (t) by a series of 3 exponential contributions where the "effective" remaining mass m 3 , resumming any contributions but the ground state and the first excited state, can be different for every correlator. Decay constants of J/ψ and ψ(2S) are proportional to Z 0 1 / √ m V1 and Z 0 2 / √ m V2 (the local-local vector 2-pt correlator corresponds to i = j = 0) and we find the hierarchy Z 0 2 / √ m V2 Z 0 1 / √ m V1 in our sets of data. There is no hope to get in the continuum limit f ψ(2S) /f J/ψ < 1 using this procedure either. Neglecting disconnected diagrams might be a source of the problem; including more flavours in the sea might help to reduce the spin breaking effects, as it was observed on quantities like f D * s /f Ds [34][35][36][37]. Besides, the width Γ(ψ(2S) → e + e − ) is smaller than its ground state counterpart Γ(J/ψ → e + e − ), that is 2.34 keV versus 5.56 keV [33]. Written in terms of the decay constant f ψ(2S) and the mass m ψ(2S) , this is a serious clue that f ψ(2S) /f J/ψ < 1 1 . Moreover, a lattice study, performed in the framework of NRQCD, gives f Υ /f Υ < 1 in the bottomium sector [38].

III. CONCLUSION
In that paper we have reported a N f = 2 lattice QCD study about the physics of quarkonia. The decay constants f ηc and f J/ψ are in the same ballpark as the two previous lattice estimations available so far in the literature, with the good news that cut-off effects seem to be limited to 10%. The issues with radial excitations are difficult to circumvent. As a first solution, the basis of operators in the GEVP analysis could be enlarged by including interpolating fields with covariant derivatives or operators of the π 2 and ρ 2 kind. But both of them suffer either from big statistical fluctuations, because of numerical cancellation among various contributions, or from the more serious conceptual problem that, in GEVP, mixing Teven and T-odd operators has no real sense. Second, m ηc(2S) /m ηc and m ψ(2S) /m J/ψ are not significantly affected by cut-off effects, so that they are 5% larger than the experimental ratios. The information that f ηc(2s) /f ηc ∼ 0.8 is interesting, as the decay constants f ηc and f ηc(2S) are hadronic inputs that govern the transitions η c → l + l − , h → η c l + l − , η c (2S) → l + l − and h → η c (2S)l + l − with a light CP-odd Higgs boson as an intermediate state 2 . Unfortunately, our result for f ψ(2S) /f J/ψ > 1 makes the picture less bright, unless one admits that there are very large spin breaking effects. Further investigation to address this issue is undoubtedly required. Nevertheless, the next step into the measurement of f η b , particularly relevant in models with a light CPodd Higgs, is underway using step scaling in masses in order to extrapolate the results to the bottom region.