(cid:2) c (cid:3) c π coupling and (cid:3) c → (cid:2) c π decay in lattice QCD

We evaluate the (cid:2) c (cid:3) c π coupling constant ( G (cid:2) c (cid:3) c π ) and the width of the strong decay (cid:3) c → (cid:2) c π in 2 + 1 ﬂavor lattice QCD on four different ensembles with pion masses ranging from 700 MeV to 300 MeV. We ﬁnd G (cid:2) c (cid:3) c π = 18 . 332 ( 1 . 476 ) stat . ( 2 . 171 ) syst . and the decay width (cid:5)((cid:3) c → (cid:2) c π ) = 1 . 65 ( 28 ) stat . ( 30 ) syst . MeV on the physical quark-mass point, which is in agreement with the recent experimental determination. the


Introduction
We have seen an immense progress on the physics of charmed baryons in the last decade and all the ground-state singlecharmed baryons and several excited states, as predicted by the quark model, have been experimentally measured [ The strong decay c → c π has been studied in Heavy Hadron Chiral Perturbation Theory [9][10][11], Light-front Quark Model [12], Relativistic Quark Model [13], nonrelativistic Quark Model [14,15], 3 P 0 Model [16] and QCD Sum Rules [17]. Most recently, Belle Collaboration has measured the decay width of c (2455) ++ as = 1.84 ± 0.04 +0.07 −0.20 MeV and that of c (2455) 0 as = 1.76 ± 0.04 +0.09 −0.21 MeV [18]. We have recently extracted the electromagnetic form factors of baryons in lattice QCD [19][20][21]. Motivated by the recent experi-mental measurements, in this work we broaden our program to include pion couplings of baryons. As a first step we evaluate the strong coupling constant c c π and the width of the strong decay c → c π in 2 + 1 flavor lattice QCD. Our aim is to utilize this calculation as a benchmark for future calculations. This work is reminiscent of Refs. [22,23] where pion-octet-baryon coupling constants have been calculated in lattice QCD.
Our work is organized as follows: In Section 2 we present the theoretical formalism of our calculations of the form factors together with the lattice techniques we have employed to extract them. In Section 3 we present and discuss our numerical results. Section 4 contains a summary of our findings.

Theoretical formulation and lattice simulations
We begin with formulating the baryon matrix elements of the pseudoscalar current, which we evaluate on the lattice to compute the pion coupling constants. The pion has a direct coupling to the axial-vector current A a where f π = 92 MeV is the pion decay constant. Taking the divergence of the axial-vector current, we find the partially conserved axial-vector current (PCAC) hypothesis where φ a is the pion field with the normalization 0|φ a (0)|π b (q) = δ ab . The matrix element of the PCAC hypothesis between baryon states yields where P a (x) =ψ(x)γ 5 τ a 2 ψ(x) is the pseudoscalar current and ψ(x) is the isospin doublet quark field. Inserting Eq. (4) into Eq.
(3), we find the baryon-baryon matrix elements of the pseudoscalar current which we use to extract G B Bπ . We use the values of pion decay constant, f π , pion mass, m π , and the quark mass, m q , on each ensemble as determined by PACS-CS [24]. The dependence on the pseudoscalar current renormalization constant cancels on the lefthand side of Eq. (5). While the matrix element in Eq. (5) is derived by a PCAC prescription we can extract the pseudoscalar matrix elements on the lattice directly by using the following ratio R(t 2 , t 1 ; p , p; ; μ) where the baryonic two-point and three-point correlation functions are respectively defined as G B PB (t 2 , t 1 ; p , p; ) with i = γ i γ 5 4 and 4 ≡ (1 + γ 4 )/2. t 1 is the time when the external pseudoscalar field interacts with a quark and t 2 is the time when the final baryon state is annihilated.
The baryon interpolating fields are chosen as where i, j, k denote the color indices and C = γ 4 γ 2 . In the large Euclidean time limit, t 2 − t 1 and t 1 a, the ratio in Eq. (6) reduces to the desired form where Q 2 = −q 2 . We measure the c c π coupling constant for both kinematical cases with B = c , B = c (denoted by G c c π ) and B = c , B = c (denoted by G c c π ).
Here we summarize our lattice setup and refer the reader to Ref. [25] for the details since we employ the same setup in this work. We have run our lattice simulations on 32 3 × 64 lattices with 2 + 1 flavors of dynamical quarks using the gauge configurations generated by the PACS-CS collaboration [24] with the nonperturbatively O(a)-improved Wilson quark action and the Iwasaki gauge action. We use the gauge configurations at β = 1.90 with the clover coefficient c S W = 1.715 having a lattice spacing of a = 0.0907(13) fm (a −1 = 2.176(31) GeV). We consider four different hopping parameters for the sea and the u, d valence quarks, κ sea , κ u,d val = 0.13700, 0.13727, 0.13754 and 0.13770, which correspond to pion masses of ∼ 700, 570, 410, and 300 MeV, respectively. We also include data with κ sea , κ u,d val = 0.13781 for mass determination.
We use the wall method which does not require to fix sink operators in advance and hence allowing us to compute all baryon channels we are interested in simultaneously. However, since the wall sink/source is a gauge-dependent object, we have to fix the gauge, which we choose to be Coulomb. We extract the baryon masses from the two-point correlator with shell source and point sink, and use the dispersion relation to calculate the energy at each momentum transfer.
Similar to our simulations in Ref. [25], we choose to employ Clover action for the charm quark. While the Clover action is subject to discretization errors of O(m q a), it has been shown that the calculations which are insensitive to a change of charm-quark mass are less severely affected by these errors [19][20][21]25,26]. Note that the Clover action we are employing here is a special case of the Fermilab heavy-quark action with c S W = c E = c B [27]. We determine the hopping parameter of the charm quark nonperturbatively as κ c = 0.1246 by tuning the spin-averaged static masses of charmonium and heavy-light mesons to their experimental values [20].
We employ smeared source and wall sink which are separated by 12 lattice units in the temporal direction. Light and charm quark source operators are smeared in a gauge-invariant manner with the root mean square radius of r l ∼ 0.5 fm and r c = r l /3 respectively. All the statistical errors are estimated via the jackknife analysis. In this work, we consider only the connected diagrams since the P 3 current is an isovector current and the relevant light quark disconnected diagrams vanish. We make our measurements on 100, 100, 200 and 315 configurations, respectively for each quark mass. In order to increase the statistics we take several different source points using the translational invariance along the temporal direction. We make momentum insertions in all directions and average over equivalent (positive and negative) momenta. Computations are performed using a modified version of Chroma software system [28] on CPU clusters along with QUDA [29,30] for propagator inversion on GPUs.

Results and discussion
Masses of the baryons in question are input parameters for form factor calculations. In Table 1, we give c and c masses for five light-quark hopping-parameter values corresponding to each light-quark mass we consider. We extrapolate the masses to the physical point by a HHχ PT procedure as outlined in Ref. [31]. Our results are compared to those reported by PACS-CS [32], ETMC [33], Briceno et al. [34] and Brown et al. [35] and to the experimental values [1] in Table 1.
Our previous determinations of the charmed-baryon masses [19,20]  In this work, we use the physical-point results to estimate systematic errors due to employing Clover action for charm  (14) quarks and our κ c tuning rather than a detailed spectroscopic analysis. Since the baryon masses only appear in kinematical terms in form factor calculations, the sensitivity of the final results to mass deviations are negligible. Considering the ∼ 2% discrepancy between our and experimental mass values, we expect that any systematic error due to charm quark would have a similar or less effect on the form factor values for which the statistical errors are much larger. We make our analysis by considering two different kinematic cases where we choose the source particle as a c or a c particle. The first case corresponds to the c → c π transition where the particle at sink, that is c , is at rest since its momentum is projected to zero due to wall smearing. The second case is the c π → c transition where c is located at the sink point.
A common practice to extract the form factors is to identify the regions where the ratio in Eq. (6) remains constant, namely forms a plateau with respect to the current-insertion time, t 1 . However, due to a finite source-sink separation, it might not always be possible to identify a clean plateau signal and an asymmetric (Gaussian smeared) source-(wall smeared) sink pair, as employed here would further affect the signal since different smearing procedures are known to cause different ground-state approaches. An ill-defined plateau range would be prone to excited state contamination which would introduce an uncontrolled systematic error. In order to check that our plateau analysis yields reliable results we compare the form factor values extracted by the plateau method to the ones extracted by a phenomenological form given as, where the first term is the form factor value we wish to extract and the coefficients b 1 , b 2 and the mass gaps 1 , 2 are regarded as free parameters. Fig. 1 shows the ratio in Eq. (6) as a function of currentinsertion time t 1 with 12a (∼ 1.09 fm) separation between the source and the sink on the heaviest quark ensemble (κ u,d = 0.13700) and for various momentum insertions. We compare the two form-factor values as extracted by the plateau method and by the phenomenological form fits. Apparent discrepancy between different fit procedures in the c π → c kinematical case hints that either the data set is unreliable or the analysis suffers from excited-state contaminations. On the other hand, the c → c π case exhibits a good agreement between a plateau and a phenomenological approach. We observe a similar behavior on the other ensembles also as shown in the Fig. 5. We utilize the phenomenological form as a cross check rather than the actual fit procedure since regression analysis has a tendency to become unstable with increased number of free parameters. As long as the plateau fit results agree with that of the phenomenological form fit's we deem the data as reliable, less prone to excited state contamination and thus trust the identified plateaux and adopt its values for form factors.
As a further check of possible excited-state contaminations, we repeat the simulations on the κ u,d = 0.13700 ensemble with a larger source-sink separation of 14 lattice units (∼ 1.27 fm). Fig. 2 shows the ratio in Eq. (6) as a function of current-insertion time for various momentum insertions with t 2 = 12 and t 2 = 14. In the case of c π → c there is a large discrepancy between the R(t 2 , t 1 ; p , p; ; μ) values of two different source-sink separations and furthermore data are systematically smaller unlike the phenomenological form fit results. This inconsistency implies that not only the c π → c case has significant excited state contamination but also the plateau and phenomenological-form fit analyses of the 12a data is unreliable. On the other hand, the 12a and 14a behavior of the c → c π case is similar and consistent with the 12a phenomenological form analysis leading us to infer that c → c π is less affected by excited-state contaminations.   13770. While all form factors have a tendency to decrease as momentum transfer increases, there is a visible correlation amongst the data corresponding to first three and second three Q 2 values. Note that a similar behavior also appears in the previous works on pseudoscalar-baryon coupling constants [22,23]. One possible source of this clustering with respect to momenta is the uncontrolled systematic errors such as discretization errors, which can be mitigated by use of finer lattices. In order to circumvent this problem one can analyze the on-axis (all momenta carried on a single axis; i.e. (p x , p y , p z ) = (0, 0, 1), (0, 0, 2) and (0, 0, 3)) data only and perform a functional-form fit to extract the values at Q 2 = 0. Such an analysis however discards useful low-momentum data which is crucial to constrain the fits. We note that although we do not rely on this method, except in the κ u,d = 0.13770 case where the signal deteriorates heavily, our results given below differ by less than 3% from those of an on-axis analysis. One other source for the clustering of data might be Lorentz symmetry breaking and hyper-cubic effects. Hyper-cubic lattice artefacts can be identified from observables extracted at a given p 2 value with different momentum combinations, e.g., the form factor evaluated with (p x , p y , p z ) = (0, 0, 3) and (2, 2, 1). We have made this test by measuring G c c π with (p x , p y , p z ) = (0, 0, 3) and (2, 2, 1) momentum combinations and on 100 configurations with κ u,d =  [36][37][38], however, we need more data with similar momentum combinations to make a conclusive analysis. Note that when the data with momentum combination (p x , p y , p z ) = (2, 2, 1) instead of that with (p x , p y , p z ) = (0, 0, 3) (or their average) is used in the Q 2 fit, the fitted results of form factors at Q 2 = 0 are only slightly affected.
We perform fits to Q 2 using pole-form ansätze, viz. a monopole form and a dipole form as given below, where the is a free pole-mass parameter. We require the extrapolated values to Q 2 = 0 using two ansätze to be as close to each other as possible since the coupling constant value should be independent of the ansatz that's used to describe the form factors. We observe that such a condition is best realized in the c → c π case.
In order to make the final consideration to quantify the systematic errors arising due to the excited-state contamination, we visit the comparison of two cases with source-sink separation values once again and compare the extrapolated coupling constants.
We show the plots of form factors with t 2 = 12a and t 2 = 14a in Fig. 3 where each data point is extracted by a plateau analysis.
We focus particularly on the c → c π case for which the ex- One important observation from the c → c π kinematical case in Fig. 3 is that the correlation amongst the data mentioned above seems to vanish when the source sink separation is increased. However, any apparent correlation might be hidden by the increased statistical uncertainty. We have performed the t 2 = 12a and t 2 = 14a analysis with the same number of ensembles and the statistical errors increase roughly by 50%. It would require at least twice as many measurements to reach a similar precision of t 2 = 12a case. Although plausible for the κ u,d = 0.13700 case, this would not be possible for lighter quark-mass ensembles since the number of gauge configurations available is limited.
Our conclusion from the above analysis is that the c → c π kinematical case with t 2 = 12a source-sink separation is less prone to excited-state contaminations and therefore we give our final results considering the c → c π kinematical case only. We will assign a systematic error of minimum 6% to the weighted averages of the coupling constants and propagate that error to the decay width in addition to the statistical errors.
We have tabulated the coupling constants as extracted on each ensemble with different functional forms in Table 2. In Fig. 4 we show the m 2 π dependence of the G c c π (Q 2 = 0). We regard the deviation arising from different ansätze used as a source of systematic error in our calculation and estimate the error by comparing the weighted average of monopole and dipole fit results to the dipole fit result on the physical point. Lower panel of Table 2 gives the results of the extrapolations to the physical point by a constant, by a linear and by a more general quadratic form in m 2 π . There is a reasonable agreement between the results of different extrapolation forms to the physical point. The weighted averages, reported on the final column of Table 2, agree well with each other.
The final value we quote for the coupling constant is, where the first error is statistical and the second one is the combined systematical error due to weighted average and excited state contamination.
Using the physical values of the baryon masses reported by the PDG [1], we evaluate the decay width given in Eq. (15) as c = 1.65 ± 0.28 stat. ± 0.30 syst. MeV, which is in agreement with the recent experimental decay width determination of different isospin states as  [18]. For comparison, we compile other theoretical determinations of the decay widths in the literature in Table 3. In general other theoretical works tend to overestimate the c decay width as compared to experiment and our lattice result.

Conclusion
In summary, we have evaluated the c c π coupling constant and the width of the strong decay c → c π in 2 + 1 flavor lattice QCD on four different ensembles with pion masses ranging from ∼ 700 to 300 MeV. A systematic analysis of different kinematical cases and the excited state contributions is given. Incorporating our results into the strong c → c π decay, we have obtained the decay width of c as ( c → c π) = 1.65(28)(30) MeV, which is in agreement with the experimental determination. Table 3 Comparison of our result with those from experiment [18], Heavy Hadron Chiral Perturbation Theory (HHχ PT) [11,10], Light-front Quark Model (LFQM) [12], Relativistic Quark Model (RQM) [13], Non-Relativistic Quark Model (NRQM) [14,15], 3 P 0 Model [16] and QCD Sum Rules (QCDSR) [17] for the decay width of c . We quote either ( ++ c → c π + ), ( 0 c → c π − ) or the isospin average. All values are given in MeV.