Decay constants of pseudoscalar mesons D,   $ D_s $  ,   $ D_s $   and vector mesons   $ D_s $  ,   $ D_s $  ,   $ D_s $   are determined from the   $ D_s $   lattice QCD at a lattice spacing   $ D_s $   fm. For vector mesons, the decay constants defined by tensor currents are given in the   $ D_s $   scheme at   $ D_s $   GeV. The calculation is performed on domain wall fermion configurations generated by the RBC-UKQCD collaborations and the overlap fermion action is used for the valence quarks. Comparing the current results with our previous results at a coarser lattice spacing   $ D_s $   fm provides a better understanding of the discretization error. We obtain   $ D_s $   with a better precision than our previous result. Combining our   $ D_s $   MeV with the total width of   $ D_s $   determined in a recent study gives a branching fraction   $ D_s $   for   $ D_s $   leptonic decay.

China

Dong-Hao Li

Ying Chen

Ming Gong

United States

Keh-Fei Liu

Zhaofeng Liu

Ting-Xiao Wang

(χQCD Collaboration) 

Charm physics with overlap fermions on 2+1-flavor domain wall fermion configurations  

We calculate the form factor   $M(q^2)$   for the Dalitz decay   $M(q^2)$   with   $M(q^2)$   being the SU(  $M(q^2)$  ) flavor singlet pseudoscalar meson. The difference among the partial widths   $M(q^2)$   at different   $M(q^2)$   can be attributed in part to the   $M(q^2)$   and quark mass dependences induced by the   $M(q^2)$   anomaly dominance.   $M(q^2)$   in both   $M(q^2)$   is well described by the single pole model   $M(q^2)$  . Combined with the known experimental results of the Dalitz decays   $M(q^2)$  , the pseudoscalar mass   $M(q^2)$   dependence of the pole parameter   $M(q^2)$   is approximated by   $M(q^2)$   with   $M(q^2)$   and   $M(q^2)$  . These results provide inputs for future theoretical and experimental studies on the Dalitz decays   $M(q^2)$  .

Chunjiang Shi

Xiangyu Jiang

Wei Sun

Form factor for Dalitz decays from / To Light Pseudoscalars  

By extracting the transition amplitudes, we give the first lattice QCD prediction of the two-body decay partial widths of the <math><msup><mn>1</mn><mrow><mo>−</mo><mo>+</mo></mrow></msup></math> charmoniumlike hybrid <math><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></math>. Given the calculated mass value <math><msub><mi>m</mi><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></msub><mo>=</mo><mn>4.329</mn><mo>(</mo><mn>36</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>GeV</mi></math>, the <math><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></math> decay is dominated by the open charm modes <math><msub><mi>D</mi><mn>1</mn></msub><mover><mi>D</mi><mo>¯</mo></mover></math>, <math><msup><mi>D</mi><mo>*</mo></msup><mover><mrow><mi>D</mi><mo>,</mo></mrow><mo>¯</mo></mover></math> and <math><msup><mi>D</mi><mo>*</mo></msup><msup><mover><mi>D</mi><mo>¯</mo></mover><mo>*</mo></msup></math> with partial widths of 258(133), 88(18), and 150(118) MeV, respectively. The coupling of <math><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></math> to <math><msub><mi>χ</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></math> plus a flavor singlet pseudoscalar is not small, but <math><msub><mi>χ</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub><mi>η</mi></math> decay is suppressed by the small <math><mi>η</mi><mo>−</mo><msup><mi>η</mi><mo>′</mo></msup></math> mixing angle. The partial width of <math><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub><mo>→</mo><msub><mi>η</mi><mi>c</mi></msub><msup><mi>η</mi><mo>′</mo></msup></math> is estimated to be around 1 MeV. We suggest experiments to search for <math><msub><mi>η</mi><mrow><mi>c</mi><mn>1</mn></mrow></msub></math> in the <math><mi>P</mi></math>-wave <math><msup><mi>D</mi><mo>*</mo></msup><mover><mi>D</mi><mo>¯</mo></mover></math> and <math><msup><mi>D</mi><mo>*</mo></msup><msup><mover><mi>D</mi><mo>¯</mo></mover><mo>*</mo></msup></math> systems. Especially, the polarization of <math><msup><mi>D</mi><mo>*</mo></msup><msup><mover><mi>D</mi><mo>¯</mo></mover><mo>*</mo></msup></math> can be used to distinguish the <math><msup><mn>1</mn><mrow><mo>−</mo><mo>+</mo></mrow></msup></math> product (total spin <math><mi>S</mi><mo>=</mo><mn>1</mn></math>) from <math><msup><mn>1</mn><mrow><mo>−</mo><mo>−</mo></mrow></msup></math> products (<math><mi>S</mi><mo>=</mo><mn>0</mn></math>).

Decays of <math><msup><mn>1</mn><mrow><mo>−</mo><mo>+</mo></mrow></msup></math> charmoniumlike hybrid using lattice QCD

We present the light-flavor quark masses and low-energy constants using the <math><mrow><mn>2</mn><mo>+</mo><mn>1</mn></mrow></math> flavor full-QCD ensembles with stout smeared-clover fermion action and Symanzik gauge action. Both the fermion and gauge actions are tadpole improved self-consistently. The simulations are performed on 11 ensembles at three lattice spacings <math><mi>a</mi><mo>∈</mo><mo>[</mo><mn>0.05</mn><mo>,</mo><mn>0.11</mn><mo>]</mo><mtext> </mtext><mtext> </mtext><mi>fm</mi></math>, four spatial sizes <math><mi>L</mi><mo>∈</mo><mo>[</mo><mn>2.5</mn><mo>,</mo><mn>5.1</mn><mo>]</mo><mtext> </mtext><mtext> </mtext><mi>fm</mi></math>, seven pion masses <math><msub><mi>m</mi><mi>π</mi></msub><mo>∈</mo><mo>[</mo><mn>135</mn><mo>,</mo><mn>350</mn><mo>]</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>, and several values of the strange quark mass. The quark mass is defined through the partially conserved axial current relation and renormalized to <math><mrow><mover><mrow><mi>MS</mi></mrow><mrow><mo>¯</mo></mrow></mover><mo>(</mo><mn>2</mn><mtext> </mtext><mi>GeV</mi><mo>)</mo></mrow></math> through the intermediate regularization independent momentum subtraction scheme. The systematic uncertainty of using the symmetric momentum subtraction scheme is also included. Eventually, we predict <math><msub><mi>m</mi><mi>u</mi></msub><mo>=</mo><mn>2.45</mn><mo>(</mo><mn>22</mn><mo>)</mo><mo>(</mo><mn>20</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>, <math><msub><mi>m</mi><mi>d</mi></msub><mo>=</mo><mn>4.74</mn><mo>(</mo><mn>11</mn><mo>)</mo><mo>(</mo><mn>09</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>, and <math><msub><mi>m</mi><mi>s</mi></msub><mo>=</mo><mn>98.8</mn><mo>(</mo><mn>2.9</mn><mo>)</mo><mo>(</mo><mn>4.7</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> with the systematic uncertainties from lattice spacing determination, continuum extrapolation and renormalization constant included. We also obtain the chiral condensate <math><msup><mi>Σ</mi><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup><mo>=</mo><mn>268.6</mn><mo>(</mo><mn>3.6</mn><mo>)</mo><mo>(</mo><mn>0.7</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> and the pion decay constant <math><mi>F</mi><mo>=</mo><mn>86.6</mn><mo>(</mo><mn>7</mn><mo>)</mo><mo>(</mo><mn>1.4</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> in the <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>2</mn></math> chiral limit, and the next-to-leading order low-energy constants <math><msub><mo>ℓ</mo><mn>3</mn></msub><mo>=</mo><mn>2.43</mn><mo>(</mo><mn>54</mn><mo>)</mo><mo>(</mo><mn>05</mn><mo>)</mo></math> and <math><msub><mo>ℓ</mo><mn>4</mn></msub><mo>=</mo><mn>4.322</mn><mo>(</mo><mn>75</mn><mo>)</mo><mo>(</mo><mn>96</mn><mo>)</mo></math>.

Zhi-Cheng Hu

Bo-Lun Hu

Ji-Hao Wang

Guoming Liu

Liuming Liu

Peng Sun

Wei Wang

Yi-Bo Yang

Dian-Jun Zhao

Quark masses and low-energy constants in the continuum from the tadpole-improved clover ensembles

On configurations with <math><mrow><mn>2</mn><mo>+</mo><mn>1</mn></mrow></math>-flavor dynamical domain-wall fermions, we calculate the RI/(S)MOM renormalization constants (RC) of overlap quark bilinears. Hypercubic (HYP) smearing is used to construct the overlap Dirac operator. We investigate the possible effects of the smearing on discretization errors in the RCs by varying the level of smearing from 0 to 1 and 2. The lattice is of size <math><msup><mn>32</mn><mn>3</mn></msup><mo>×</mo><mn>64</mn></math> and with lattice spacing <math><mrow><mn>1</mn><mo>/</mo><mi>a</mi><mo>=</mo><mn>2.383</mn><mo>(</mo><mn>9</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>GeV</mi></mrow></math>. The RCs in the <math><mover><mi>MS</mi><mo>¯</mo></mover></math> scheme at 2 GeV are given at the end, with the uncertainty of <math><msub><mi>Z</mi><mi>T</mi></msub></math> reaching <math><mo>≤</mo><mn>1</mn><mo>%</mo></math> for the tensor current. Results of the renormalized quark masses and hadron matrix elements show that the renormalization procedure suppresses the <math><mo>∼</mo><mn>30</mn><mo>%</mo></math> difference of the bare quantities with or without HYP smearing into the 3%–5% level.

Yujiang Bi

Fangcheng He

RI/(S)MOM renormalizations of overlap quark bilinears with different levels of hypercubic smearing

We perform the first lattice study on the mixing of the isoscalar pseudoscalar meson <math><mi>η</mi></math> and the pseudoscalar glueball <math><mi>G</mi></math> in the <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>2</mn></math> QCD at the pion mass <math><msub><mi>m</mi><mi>π</mi></msub><mo>≈</mo><mn>350</mn><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>. The <math><mi>η</mi></math> mass is determined to be <math><msub><mi>m</mi><mi>η</mi></msub><mo>=</mo><mn>714</mn><mo>(</mo><mn>6</mn><mo>)</mo><mo>(</mo><mn>16</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>. Through the Witten-Veneziano relation, this value can be matched to a mass value of <math><mo>∼</mo><mn>981</mn><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> for the SU(3) counterpart of <math><mi>η</mi></math>. Based on a large gauge ensemble, the <math><mi>η</mi><mo>−</mo><mi>G</mi></math> mixing energy and the mixing angle are determined to be <math><mo>|</mo><mi>x</mi><mo>|</mo><mo>=</mo><mn>107</mn><mo>(</mo><mn>15</mn><mo>)</mo><mo>(</mo><mn>2</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> and <math><mrow><mo>|</mo><mi>θ</mi><mo>|</mo><mo>=</mo><mn>3.46</mn><mo>(</mo><mn>46</mn><mo>)</mo><mo>°</mo></mrow></math> from the <math><mi>η</mi><mo>−</mo><mi>G</mi></math> correlators that are calculated using the distillation method. We conclude that the <math><mi>η</mi><mo>−</mo><mi>G</mi></math> mixing is tiny and the topology induced interaction contributes most of <math><mi>η</mi></math> mass owing to the QCD <math><msub><mi>U</mi><mi>A</mi></msub><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></math> anomaly.

Feiyu Chen

Renqiang Zhang

<math><mi>η</mi></math>-glueball mixing from <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>2</mn></math> lattice QCD

We present the first theoretical prediction of the partial decay width of the process <math><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><mi>γ</mi><msub><mi>η</mi><mn>1</mn></msub></math>, where <math><msub><mi>η</mi><mn>1</mn></msub></math> is the lightest flavor singlet <math><msup><mn>1</mn><mrow><mo>−</mo><mo>+</mo></mrow></msup></math> hybrid meson. Our <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>2</mn></math> lattice QCD calculation at <math><msub><mi>m</mi><mi>π</mi></msub><mo>≈</mo><mn>350</mn><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math> results in the <math><msub><mi>η</mi><mn>1</mn></msub></math> mass <math><msub><mi>m</mi><msub><mi>η</mi><mn>1</mn></msub></msub><mo>=</mo><mn>2.23</mn><mo>(</mo><mn>4</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>GeV</mi></math> and the related electromagnetic form factors <math><msub><mi>M</mi><mn>1</mn></msub><mo>(</mo><mn>0</mn><mo>)</mo><mo>=</mo><mo>−</mo><mn>4.73</mn><mo>(</mo><mn>74</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>, <math><msub><mi>E</mi><mn>2</mn></msub><mo>(</mo><mn>0</mn><mo>)</mo><mo>=</mo><mn>1.18</mn><mo>(</mo><mn>22</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>, which give <math><mi>Γ</mi><mo>(</mo><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><mi>γ</mi><msub><mi>η</mi><mn>1</mn></msub><mo>)</mo><mo>=</mo><mn>2.04</mn><mo>(</mo><mn>61</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>eV</mi></math>. These form factors can be applied to the physical <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>3</mn></math> case, where there should be two hybrid mass eigenstates <math><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mi>l</mi><mo>)</mo></mrow></msubsup></math> and <math><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mi>h</mi><mo>)</mo></mrow></msubsup></math> due to the singlet-octet mixing. It is shown that the ratio of the branching fractions <math><mi>Br</mi><mo>(</mo><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><mi>γ</mi><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mi>l</mi><mo>,</mo><mi>h</mi><mo>)</mo></mrow></msubsup><mo>→</mo><mi>γ</mi><mi>η</mi><msup><mi>η</mi><mo>′</mo></msup><mo>)</mo></math> is inversely proportional to the ratio of the total widths of <math><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mi>l</mi><mo>,</mo><mi>h</mi><mo>)</mo></mrow></msubsup></math>. Given our results and the mixing angle derived by a previous lattice study, whether <math><msub><mi>η</mi><mn>1</mn></msub><mo>(</mo><mn>1855</mn><mo>)</mo></math> is assigned to be <math><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></msubsup></math> or <math><msubsup><mi>η</mi><mn>1</mn><mrow><mo>(</mo><mi>h</mi><mo>)</mo></mrow></msubsup></math>, the observed branching fraction <math><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><msub><mi>η</mi><mn>1</mn></msub><mo>(</mo><mn>1855</mn><mo>)</mo><mo>→</mo><mi>γ</mi><mi>η</mi><msup><mi>η</mi><mo>′</mo></msup></math> implies a very large coupling of the octet <math><msub><mi>η</mi><mn>1</mn></msub></math> to <math><mi>η</mi><msup><mi>η</mi><mo>′</mo></msup></math>. This should be investigated in future studies.

<math><msup><mn>1</mn><mrow><mo>−</mo><mo>+</mo></mrow></msup></math> hybrid meson in <math><mi>J</mi><mo>/</mo><mi>ψ</mi></math> radiative decays from lattice QCD

The large radiative production rate for pseudoscalar mesons in the <math><mi>J</mi><mo>/</mo><mi>ψ</mi></math> radiative decay remains elusive. We present the first lattice QCD calculation of partial decay width of <math><mi>J</mi><mo>/</mo><mi>ψ</mi></math> radiatively decaying into <math><msub><mi>η</mi><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></msub></math>, the SU(2) flavor singlet pseudoscalar meson, which confirms QCD <math><mrow><msub><mi>U</mi><mi>A</mi></msub></mrow><mo>(</mo><mn>1</mn><mo>)</mo></math> anomaly enhancement to the coupling of gluons with flavor singlet pseudoscalar mesons. The lattice simulation is carried out using <math><msub><mi>N</mi><mi>f</mi></msub><mo>=</mo><mn>2</mn></math> lattice QCD gauge configurations at the pion mass <math><msub><mi>m</mi><mi>π</mi></msub><mo>≈</mo><mn>350</mn><mtext> </mtext><mtext> </mtext><mi>MeV</mi></math>. In particular, the distillation method has been utilized to calculate light quark loops. The results are reported here with the mass <math><mrow><msub><mrow><mi>m</mi></mrow><mrow><msub><mrow><mi>η</mi></mrow><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></msub></mrow></msub><mo>=</mo><mn>718</mn><mo>(</mo><mn>8</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>MeV</mi></mrow></math> and the decay width <math><mrow><mi>Γ</mi><mo>(</mo><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><mi>γ</mi><msub><mrow><mi>η</mi></mrow><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></msub><mo>)</mo><mo>=</mo><mn>0.385</mn><mo>(</mo><mn>45</mn><mo>)</mo><mtext> </mtext><mtext> </mtext><mi>keV</mi></mrow></math>. By assuming the dominance of <math><msub><mi>U</mi><mi>A</mi></msub><mo>(</mo><mn>1</mn><mo>)</mo></math> anomaly and flavor singlet-octet mixing angle <math><mrow><mi>θ</mi><mo>=</mo><mo>−</mo><mn>24.5</mn><mi>°</mi></mrow></math>, the production rates for the physical <math><mi>η</mi></math> and <math><msup><mi>η</mi><mo>′</mo></msup></math> in <math><mi>J</mi><mo>/</mo><mi>ψ</mi></math> radiative decay are predicted to be <math><mn>1.15</mn><mo>(</mo><mn>14</mn><mo>)</mo><mo>×</mo><msup><mn>10</mn><mrow><mo>−</mo><mn>3</mn></mrow></msup></math> and <math><mn>4.49</mn><mo>(</mo><mn>53</mn><mo>)</mo><mo>×</mo><msup><mn>10</mn><mrow><mo>−</mo><mn>3</mn></mrow></msup></math>, respectively, which agree well with the experimental measurement data. Our study manifests the potential of lattice QCD studies on the light hadron production in <math><mi>J</mi><mo>/</mo><mi>ψ</mi></math> radiative decays.

Ning Li

Radiative Decay Width of <math><mrow><mi>J</mi><mo>/</mo><mi>ψ</mi><mo>→</mo><mi>γ</mi><msub><mrow><mi>η</mi></mrow><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></msub></mrow></math> from <math><mrow><msub><mrow><mi>N</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>=</mo><mn>2</mn></mrow></math> Lattice QCD

The S-wave <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> scattering in the isospin <math><mi>I</mi><mo>=</mo><mn>0</mn><mo>,</mo><mn>1</mn></math> channels is studied in <math><msub><mrow><mi>N</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>=</mo><mn>2</mn></math> lattice QCD at <math><msub><mrow><mi>m</mi></mrow><mrow><mi>π</mi></mrow></msub><mo>≈</mo><mn>350</mn></math> MeV. It is observed that the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> interaction is repulsive in the <math><mi>I</mi><mo>=</mo><mn>1</mn></math> channel when the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> energy is near the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> threshold. In contrast, the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> interaction in the <math><mi>I</mi><mo>=</mo><mn>0</mn></math> channel is definitely attractive in a wide range of the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> energy. This is consistent with the isospin assignment <math><mi>I</mi><mo>=</mo><mn>0</mn></math> for <math><msubsup><mrow><mi>T</mi></mrow><mrow><mi>c</mi><mi>c</mi></mrow><mrow><mo>+</mo></mrow></msubsup><mo>(</mo><mn>3875</mn><mo>)</mo></math>. By analyzing the components of the <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> correlation functions, it turns out that the quark diagram responsible for the different properties of <math><mi>I</mi><mo>=</mo><mn>0</mn><mo>,</mo><mn>1</mn></math> <math><mi>D</mi><msup><mrow><mi>D</mi></mrow><mrow><mo>⁎</mo></mrow></msup></math> interactions can be understood as the charged ρ meson exchange effect. This observation provides direct information on the internal dynamics of <math><msubsup><mrow><mi>T</mi></mrow><mrow><mi>c</mi><mi>c</mi></mrow><mrow><mo>+</mo></mrow></msubsup><mo>(</mo><mn>3875</mn><mo>)</mo></math>.

Siyang Chen

Relevant Scattering From = 2 lattice QCD

In this work, we generate gauge configurations with   $ N_f = 2 $   dynamical charm quarks on anisotropic lattices. The mass shift of   $ N_f = 2 $   and   $ N_f = 2 $   charmonia due to the charm quark annihilation effect can be investigated directly in a manner of unitary theory. The distillation method is adopted to treat the charm quark annihilation diagrams at a very precise level. For   $ N_f = 2 $   charmonia, the charm quark annihilation effect barely changes the   $ N_f = 2 $   mass, but lifts the   $ N_f = 2 $   mass by approximately 3–4 MeV. For   $ N_f = 2 $   charmonia, this effect results in positive mass shifts of approximately 1 MeV for   $ N_f = 2 $   and   $ N_f = 2 $  , but decreases the   $ N_f = 2 $   mass by approximately 3 MeV. We did not obtain a reliable result for the mass shift of   $ N_f = 2 $  . In addition, we observed that the spin averaged mass of the spin-triplet   $ N_f = 2 $   charmonia is in good agreement with   $ N_f = 2 $  , as expected by the non-relativistic quark model and measured by experiments.