Reply to the Referees – CPC-2021-0466

In the distillation smearing, the number of eigenvector is 50. How this number is chosen?

The number of eigenvectors $N$ is chosen based on two following facts. First, it is a compromise between the numerical expenses and enough physical information of interest. The larger the value of $N$, the more complete the information of the short range physics from the operators is kept, but more eigenvectors mean more solutions of the linear equations of the fermion matrix. It is usually conjectured that the higher excited states are more related to the short range physics. Since we are focusing on the lowest states of charmonia, we do not include too many eigenvectors. In the mean time, $N$ should not be too small or the physical information relevant to lower states would be lost. Secondly, $N$ has a direct connection with the size of the smeared quark fields, that is to say, the larger $N$ means a narrower smearing size. In order for the charmonium interpolating fields to couple most with the lowest several states, $N$ should not be too large. As shown in Fig.1 in the manuscript, we have tried various values of $N=10,30,50,100$, the profile functions $\Psi(r)$ at $N=10$ and $N=30$ do not approach to zero at $r\to La_{s}/2$, which implies the smearing sizes are too large. We tried to calculate the charmonium spectrum with $N=100$ and $N=50$ using a small portion of the gauge ensemble and found that the extracted ground and first excited state masses are compatible with each other. Therefore, we choose $N=50$ in practice. It should be noted that the usual choice of $N$ is around 60 for light quarks in the literatures. Since the spatial size of charmonia is expected to be smaller than hadrons made up of light quarks, we think $N=50$ is safe for charmonia.

It was mentioned several times that mixing with a lower mass glueball will result in positive mass shift in the charmonium mass. This is counter-intuitive. Shouldn’t the mixing with a lower mass guleball result in a lower charmonium mass and thus the mass shift due to the annihilation diagram be negative ?

For a two state mixing model, the Hamiltonian can be written as

using the basis of the unmixed two states. Let $\lambda_{1,2}$ be the two eigenvalues of $H$ with $\lambda_{1}>\lambda_{2}$, the secular equation $(E_{1}-\lambda_{i})(E_{2}-\lambda_{i})=x^{2}>0$ means $\lambda_{1}>\mathrm{max}\{E_{1},E_{2}\}$ and $\lambda_{2}<\mathrm{min}\{E_{1},E_{2}\}$. That is to say, the mixing will push the energy of the higher state upward and push that of the lower state downward. If a charmonium can mix with a lighter glueball, the mass of the charmonium should be lifted a little bit.

It was argued that for the scalar and psudoscalar channels there are possible contributions from glueballs, therefore an extra term is added as in Eq.(32). As mentioned by the authors, there is also possible glueball mixing in the 2++ channel. Has Eq.(32) been tried for 2++ ?

We tried Eq.(32) for 2++ and found consistent result. Since the data of the ratio of 2++ can already be well described without the extra term in our current fitting range, adding the exponential term only increases the uncertainties of the fitting parameters.

Some discussion about the effects of missing u, d, s quarks should be enlightening.

We have added some discussions at the end of the second paragraph in the summary part of the revised manuscript.

In this paper, the mass shifts are obtained from the ratio of the disconnected part to the connected part. Is it possible to calculate the mass from the full correlation function and the connected part respectively, and then compare them to get the mass shift?

Actually, we tried both methods.

Due to the inclusion of contribution of the charm annihilation diagrams, the full correlation functions become much more noisy in the time range where the ground state has not dominated (monitored by the effective mass plots). In the pseudoscalar channel, we observed that, in the early time range, the effective mass of the full correlation functions is a little higher than that of the correlation functions with only the connected part. However it is improper to take this difference to be exactly the mass difference of the ground states caused by the disconnected parts, since the ground state cannot be disentangled reliably from the contamination from higher states in the early time range. For other channels, the differences of the effective masses of full and connected correlation function are almost non-detectable within errors.

The advantage of using the ratios to derive the mass shifts is that the ratio of the full correlation function to the connected part can suppress largely the higher state contamination. This is especially useful when the mass difference is very tiny in comparison with the absolute mass values.

page 1, line 46, right column, ”light sea quarks is ignored”, does it mean the light u,d,s quarks are all ignored, or only the u,d are ignored? It confuses me.

The light u,d and s quarks are ignored. We write it more clearly in the revised manuscript.

In fig.2, a linear fitting is used for $1^{--},1^{++},2^{++},1^{+-}$ channels, and $0^{-+}$, $0^{++}$ with an exponential term added. However, there are two different linear ranges for $1^{--}$ channel, like [4,9] , [8,12] or [8,14]. Finally, only the small region is considered, why? As I see in $1^{++}$ channel, there exits an evident linear region in [8,14]. In a word, the second linear behaviour in 1- - channel puzzles me.

The large errors of the ratio function $R(t)$ come from the disconnected diagram which becomes very noisy beyond $t/a_{t}\sim 10$. In the $1^{--}$ channel, $R(t)$ shows up good linear behavior in the time interval $[4,9]$ (very close to a constant) on both ensemble E1 and E2, so we attribute the behavior of $R(t)$ in the interval $[9,14]$ to be mostly the statistical fluctuations. Actually, the data points in this interval deviate from the linear behaviors determined from the interval $[4,9]$ by less than $2\sigma$. If we include the data points in the interval [10,14] to the fit, the central value does not change much with a mildly larger $\chi^{2}/\mathrm{dof}$. Actually on E2, all the data points in the interval $[4,14]$ converge to a single line. Since the fitted $\delta m$ in $1^{--}$ channel is very small, we do not quote the value in our conclusion but only say that the charm annihilation effects in $1^{--}$ are negligible.

We add a relevant discussion on this in the revised manuscript (in the second paragraph under Eq.(31)).

In fig.2 and fig.3, more clearly in Table.III, the two ensembles have the same lattice spacings, then it is natural to fit the values in the same time slice. In particularly, for 1- - channel in Table.III, a region [4,9] is chosen for E1, but [4,14] for E2. I am not sure whether it leads to the abnormal errors in Eq.(33), where the error in 0-+(E1) is far smaller than 0-+(E2), and error of 1+-(E1) is three times bigger than the 1+-(E2).

In our data analysis procedure, we fixed the upper bound ($t_{\mathrm{max}}$) of the fitting window and let the lower bound ($t_{\mathrm{min}}$) vary. In the $t_{\mathrm{min}}$ range where the central value of $\delta m_{1}$ is stable, we presented in the manuscript the result with the smallest $\chi^{2}/\mathrm{dof}$. While the central value of $\delta m_{1}$ is stable, its error decreases when $t_{\mathrm{min}}$ goes lower.

We have made an explanation to the different fitting ranges for the case of $1^{--}$ in response to the Question no.2.

In the $0^{-+}$ channel, the fit at $t_{\mathrm{min}}=7$ has the smallest $\chi^{2}/\mathrm{dof}$ on E2, therefore the error of $\delta m_{1}$ is larger than that at $t_{\mathrm{min}}=5$ on E1 (see Table III in the previous manuscript). As shown in Table III, for almost all channels, the fitting ranges are more or less the same for E1 and E2 ensembles. If we lower $t_{\mathrm{min}}$ to be 6 on E2, we get $\delta m_{1}(0^{-+})=3.1(2)$ MeV, which is compatible with 3.0(1) MeV on E1 at $t_{\mathrm{min}}=5$, but now the $\chi^{2}/\mathrm{dof}$ increases from 0.93 to 1.4. Since this $\chi^{2}/\mathrm{dof}$ is acceptable and comparable with 1.5 on E1 at $t_{\mathrm{min}}=5$, the result is updated to that of the new fit in the revised manuscript,

The abnormal errors of $\delta m_{1}(1^{+-})$ are due to our carelessness in the presentation of the previous manuscript. The original value of $\delta m_{1}(1^{+-})$ is 0.001002(73), but was presented to be 0.0010(7) by mistake in Table III (and accordingly 1.0(7) MeV in Eq.(33)) of the previous manuscript.

We have corrected them in the revised manuscript. We thank the Reviewer for bringing this to our attention.

For the deviation in Eq.(46), I suspect that the 2P states are not extracted exactly ,rather than precisely. To extract a higher excited-state, ones usually need to adapt a variety of different operators which are mixed with different excited-states. Though, it is feasible to choose the same gamma operators but with different smearing procedures, it is usually highly correlated with each other, as seen in fig.4. Three operators with r/a=0,3,6, may not enough to extract the 2P states completely. This is just a personal opinion.

The Reviewer’s comment is reasonable. We have tried more operators with different r/a and we have also tried 4-dimensional variation, but no better results were obtained. Our construction of the operators are relatively simple and these operators are correlated. This should be improved in the future.

The hyperfine splitting is 60 MeV, which is far smaller than PDG and also smaller than 80(2) from Hadron Spectrum Collaboration. The authors owe it to the clover terms of the fermion action and present a detailed discussion. However, the ensembles in the work have ignored the light sea quarks, Is it possible to lead to this discrepancy? The authors could present more comments on the absence of the light sea quark effects.

We have added some discussions in the summary of the revised manuscript.

In fig.6, I am confused by two $\chi^{2}/\mathrm{dof}$ for a two-states fit. Shouldn’t there be only one $\chi^{2}/\mathrm{dof}$ for such fit: c(t)=c1*(exp(-m1*t)+exp(-m1*(T-t)))+c2*(exp(-m2*t)+exp(-m2*(T-t))) ?

Sorry for the misleading. We fitted the two correlation functions for the ground state and the first excited state respectively, thus we have two $\chi^{2}/\mathrm{dof}$’s. For each correlation function we use a two-state fit. We have updated these figures in the revised manuscript to make them more readable.

The deviation of mass shift of $1^{++}$ in two ensembles are more than 3 sigma. Only difference of the ensembles is the bare quark mass. The consistencies in other channels seem to indicate the annihilation effect is not sensitive to the bare quark mass, which is not true for $1^{++}$. Are there any physical reasons and explanations for this channel?

Thank the Reviewer for pointing out this. Actually we have no idea about this discrepancy yet. In response to the Reviewer’s comment, we checked our fitting procedure by varying the fitting window. See Fig. 1 of this reply, for ensemble E1 if we fit data in the range [5,14], $\delta m_{1}(1^{++})$ is compatible with the result of ensemble E2 within $2\sigma$. Its $\chi^{2}$/dof(=1.2) is acceptable but is larger than that if we use the range [7,14]. On the other hand, there are obvious deviations from the fitting curve at large time slices for the fitting range [5,14] (see right panel of the figure). Thus we choose [7,14] to give the final result. We add a comment on this in the revised manuscript.

page 9, Eq.(43), $W_{1}$ and $W_{2}$ seem to be introduced for the first time. Its relation with $W_{in}$, defined in Eq.(40), should be identified.
$W_{1}$ and $W_{2}$ in Eq.(43) have no direct relation with $W_{in}$ and are parameters introduced to do the two-mass fit. We clarify this below Eq.(43) in the revised manuscript.