Charmonium ground and excited states at finite temperature from complex Borel sum rules

TL;DR

This paper uses complex Borel sum rules and maximum entropy methods to analyze charmonium spectral functions at finite temperature, successfully identifying both ground and excited states and their thermal modifications near the critical temperature.

Contribution

It introduces a novel approach combining complex Borel sum rules with maximum entropy to extract excited charmonia states, surpassing previous QCD sum rule limitations.

Findings

Successfully identified excited and ground charmonium states.

Observed near-simultaneous melting of all states at critical temperature.

Provided detailed spectral function modifications at finite temperature.

Abstract

Charmonium spectral functions in vector and pseudoscalar channels at finite temperature are investigated through the complex Borel sum rules and the maximum entropy method. Our approach enables us to extract the peaks corresponding to the excited charmonia, $\psi^\prime$ and $\eta_c^\prime$, as well as those of the ground states, $J/\psi$ and $\eta_c$, which has never been achieved in usual QCD sum rule analyses. We show the spectral functions in vacuum and their thermal modification around the critical temperature, which leads to the almost simultaneous melting (or peak disappearance) of the ground and excited states.

Tables

Table 1. Table 1: Values and respective uncertainties of the condensates and other parameters used for evaluating the OPE of Eq. ( 1 ).

${\overline{m}}_c\left({\overline{m}}_c\right)$	$1.273\pm 0.006$ McNeile et al. (2010)
$⟨\frac{{\alpha }_s}{\pi }{G}^2⟩$	$0.012\pm 0.0036{\mathrm{GeV}}^4$ Shifman et al. (1979b); Colangelo and Khodjamirian (2001)
${\Lambda }_{QCD}$	$0.213\pm 0.008$ Bethke (2013)

Table 2. Table 2: Positions of the charmonia peaks in vacuum, extracted from CBSR and MEM. The corresponding spectral functions are shown in Fig.  2 .

	V		PS
	1st [GeV]	2nd [GeV]	1st [GeV]	2nd [GeV]
This work	2.954	3.624	2.900	3.512
Experiment	3.096	3.686	2.983	3.639

Table 3. Table 3: Lower limits of the melting temperature determined by the MEM error estimation.

V		PS
1st peak [$T_c$]	2nd peak [$T_c$]	1st peak [$T_c$]	2nd peak [$T_c$]
1.01	0.99	1.01	0.99

Equations

$$G^J(M^2;T)$$

$$G^J(M^2;T)=\frac{1}{M^2}{\int }_0^{\infty }{\mathrm{e}}^{-s/M^2}\rho (s;T)ds.$$

$${\varphi }_b(T)$$

$${\varphi }_c(T)$$

$$⟨\frac{{\alpha }_s}{\pi }{G}^{a\mu \sigma }{G}_{\sigma }^{a\nu }{⟩}_T$$

$$G_0(T)$$

$$G_2(T)$$

$$\frac{|d^4({\mathcal{M}}^2;T=0)|}{|G^J({\mathcal{M}}^2;T=0)|}<0.1,$$

$$\begin{split}K^{\mathrm{R}}(\mathcal{M}^{2};s)&\equiv\mathrm{Re}\Bigl{[}\,\frac{1}{\mathcal{M}^{2}}\,\mathrm{e}^{-s/\mathcal{M}^{2}}\,\Bigr{]}\\ &=\frac{1}{M^{2}}\mathrm{e}^{-(\cos{\theta}/M^{2})s}\cos{\bigl{[}\,(\sin{\theta}/M^{2})s-\theta\,\bigr{]}},\end{split}$$

$$\begin{split}K^{\mathrm{I}}(\mathcal{M}^{2};s)&\equiv\mathrm{Im}\Bigl{[}\,\frac{1}{\mathcal{M}^{2}}\,\mathrm{e}^{-s/\mathcal{M}^{2}}\,\Bigr{]}\\ &=\frac{1}{M^{2}}\mathrm{e}^{-(\cos{\theta}/M^{2})s}\sin{\bigl{[}\,(\sin{\theta}/M^{2})s-\theta\,\bigr{]}}.\end{split}$$

$$\frac{\cos \theta }{M^2}>{\tilde{r}}_c.$$

Full text

I Introduction

The particle consisting of a heavy quark and an antiquark $Q\bar{Q}$, “quarkonium”, has been a suitable target to study dynamics of QCD at short distance due to its large mass Novikov et al. (1978). Furthermore, it was suggested and believed that quarkonia in an extremely hot and dense matter, quark gluon plasma (QGP), dissolve due to the color Debye screening caused by light deconfined quarks, so that such an event itself can be regarded as a signal for the existence of QGP Matsui and Satz (1986). Experimentally, this is indeed observed as the quarkonium suppression in heavy ion collisions at Relativistic Heavy Ion Collider (RHIC) at BNL and Large Hadron Collider (LHC) at CERN. On the theoretical side, this phenomenon can be understood as the thermal modification of $Q\bar{Q}$ potential Matsui and Satz (1986); Karsch et al. (1988); Karsch and Satz (1991); Digal et al. (2001); Wong (2005); Alberico et al. (2005); Satz (2006); Mócsy and Petreczky (2008, 2007); Laine et al. (2007); Rothkopf et al. (2012) or the disappearance of peaks in the spectral function by the finite temperature effects. The temperatures at which the peaks completely disappear, “melting temperatures” 111Actually, the “melting temperature” is not a well defined concept because the quarkonium wave function is gradually broadened by thermal effects, and the peak in the spectral function does not necessarily show abrupt disappearance. Therefore, in this work, we use this term just as a qualitative guideline concept., are one of the targets which should ideally be calculated from the first principles of QCD. Recently, lattice QCD simulations with the maximum entropy method (MEM) Asakawa et al. (2001) enable us to check such a spectral function deformation and to estimate the melting temperature. For instance, it was shown in such an approach that the lowest charmonium states ($J/\psi$ and $\eta_{c}$) survive above the temperature of $1.5\,T_{c}$ Asakawa and Hatsuda (2004) (cf. Umeda et al. (2005); Datta et al. (2004); Iida et al. (2006); Jakovác et al. (2007); Ding et al. (2012); Borsanyi et al. (2014); Ikeda et al. (2017)). Similarly, QCD sum rules Shifman et al. (1979a, b) can study quarkonium suppression by incorporating finite temperature effects through the QCD condensates. Initially, the ground states of charmonia and bottomonia were analyzed by assuming a specific functional form for the spectral functions Morita and Lee (2008a, b); Song et al. (2009); Dominguez et al. (2010); Morita and Lee (2010); Dominguez et al. (2011, 2013); Kim et al. (2016); Kim and Lee (2017), while recently they were reanalyzed without such assumptions by the help of the MEM Gubler et al. (2011); Suzuki et al. (2013).

In the recent heavy-ion collision experiments, the suppression of charmonium excited states ($\psi^{\prime}$ or $\psi(2S)$) was also observed Khachatryan et al. (2014); Sirunyan et al. (2017); Adam et al. (2016a, b); Adare et al. (2017); Aaij et al. (2016). The experimental results show that the yield of the excited state is more strongly suppressed when more nucleons participate in the collision. This indicates that in QGP the excited state melts at lower temperature than the ground state since the excited state suffers from the Debye screening more strongly due to the larger system size compared to the ground state. If this is the case, the second peak in the spectral function should disappear at lower temperature than the first peak also in a theoretical calculation. Actually some model studies demonstrate such a situation Matsui and Satz (1986); Karsch et al. (1988); Karsch and Satz (1991); Digal et al. (2001); Wong (2005); Alberico et al. (2005); Satz (2006); Mócsy and Petreczky (2008, 2007). Also the studies of bottomonia at finite temperature using QCD sum rules with MEM find an indirect proof indicating the same effect, without explicitly reproducing the peaks corresponding to the excited sates Suzuki et al. (2013).

In general, it is a challenging problem to extract information on excited states (namely, second or higher peaks in the spectral function) from QCD sum rules. Therefore, the thermal modification of excited states have so far not been discussed within the usual Borel type QCD sum rules, although there exist many results for the ground states. Although MEM for QCD sum rules is a powerful tool to investigate the structure of spectral function, which has already been applied to various systems Gubler and Oka (2010); Gubler et al. (2011); Ohtani et al. (2011); Suzuki et al. (2013); Ohtani et al. (2013); Gubler and Ohtani (2014); Suzuki et al. (2016); Ohtani et al. (2016), it was difficult to reproduce the second peak with statistical significance from Borel sum rules. However, by using the combination of QCD sum rules on the complex Borel plane (CBSR) and MEM, it has recently become possible to extract excited states due to the improved resolution Araki et al. (2014). In Ref. Araki et al. (2014), the excited $\phi$ meson peak in vacuum was reproduced with a mass consistent with the experimental value. The reproduction of excited charmonia from QCD sum rules is also one of main purposes of this paper.

In this work, we study charmonia in the vector and pseudoscalar channels at finite temperature by using CBSR with MEM, which has the ability to reproduce both the first and second peaks. Then we check whether there is a possible difference in the melting behaviors between the ground and excited states.

This paper is organized as follows. In Section II, we explain the CBSR for charmonia at finite temperature. In Section III, we show their spectral functions in vacuum and at finite temperature and discuss their thermal (melting) behavior. Section IV is devoted to our conclusion and outlook.

II Formalism

Let us here introduce our formalism of CBSR with MEM for the analysis of charmonium at finite temperature. Since CBSR is a generalization of real Borel sum rules even at finite temperature, we can take the same Borel sum rule formulation as done in the previous works, where we replace real Borel masses by complex ones, to construct CBSR. Here we will briefly summarize our formulation; a more detailed explanation can be found in Refs. Gubler et al. (2011); Suzuki et al. (2013). For specifying the domain of the complex Borel plane to be used in the MEM analysis, we propose an updated criterion.

II.1 The form of Borel sum rules

We calculate the correlation function for a meson system consisting of charm quarks with a mass $m$ that is larger than the typical QCD scale. Here, the dimensionless correlation functions in momentum space are defined by $\tilde{\Pi}^{V}(q^{2})\equiv\Pi_{\mu}^{V,\mu}(q)/(-3q^{2})$ and $\tilde{\Pi}^{P}\equiv\Pi^{P}(q)/q^{2}$ for vector and pseudoscalar channels, respectively. These can be calculated by the operator product expansion (OPE), and its Borel-transformed OPE with a real variable $M^{2}$ takes the following shortened form:

[TABLE]

where $\nu$ is a dimensionless parameter defined as $\nu=4m^{2}/M^{2}$. The superscript $J$ distinguishes the channel, vector or pseudoscalar. With this expression we construct sum rules at finite temperature as follows:

[TABLE]

Let us explain each part of Eq. (1). $\alpha_{s}(\nu)$ is a strong running coupling constant evaluated at the Borel mass scale $M$. The functions, $A^{J},\ a^{J},\ b^{J}\ \mbox{and}\ c^{J}$ are the Wilson coefficients. Their explicit form can be found in Ref. Morita and Lee (2010). The first and second terms come from the usual perturbative result up to the first order in $\alpha_{s}$. The $\phi(T)^{\prime}s$ are the condensates defined as

[TABLE]

where $G_{0}(T)$ and $G_{2}(T)$ are the scalar and twist-2 gluon condensates with dimension 4 at finite temperature. They are defined as the scalar and spin-2 components of the gauge independent gluon tensor:

[TABLE]

where the expression $\langle O\rangle_{T}$ means the Boltzmann average defined as $\langle O\rangle_{T}=\mathrm{Tr}(\mathrm{e}^{-H/T}O)/\mathrm{Tr}(\mathrm{e}^{-H/T})$. $u$ is the four velocity vector of the medium whose norm is unity.

II.2 Finite temperature effects

In this work we assume that all temperature dependence of the correlator enters through that of the gluon condensates. This assumption is valid as long as the temperature is smaller than the OPE separation scale, which is of the order of the Borel mass $M\sim 1,\mathrm{GeV}$. To proceed further, we therefore need to calculate the temperature dependences of the gluon condensates. In our formulation, the strategy proposed in Refs. Morita and Lee (2008a, b) is employed, in which we express the two gluon condensates as follows:

[TABLE]

where $\epsilon$, $p$ and $\alpha_{s}$ are the energy density, pressure and strong coupling constant, respectively. Thus by estimating their dependences on temperature by using quenched (pure Yang-Mills) lattice QCD simulations Boyd et al. (1996); Kaczmarek et al. (2004), we finally get the temperature dependences of the condensates. Actually the value of the critical temperature of the chiral phase transition in the quenched approximation is about $260\,\mathrm{MeV}$, while the full QCD result leads to the cross-over temperature being $145-165\,\mathrm{MeV}$ Borsanyi et al. (2010); Bazavov et al. (2012). Such a quantitative difference caused by the quenched approximation should be discussed in future works.

Let us here briefly discuss the possibility of extending our approach to full QCD, which would include dynamical quarks. Eqs. (6) and (7) in fact are derived by matching the trace part and the symmetric traceless part of the energy momentum tensor, written either in terms of thermodynamic variables or the basic degrees of freedom of QCD. While in the quenched approximation, the energy momentum tensor of QCD can be written only with gluon fields, it will have terms such as $m_{q}\langle\bar{q}q\rangle$ and $\langle\bar{q}\gamma^{\mu}D^{\nu}q\rangle$, which will hence modify Eqs. (6) and (7). Once the temperature dependences of all these operators (and $\alpha_{s}$) are known in full QCD, it will become possible to perform a QCD sum rule calculation that goes beyond the quenched approximation.

II.3 The choice of the domain in the complex Borel mass plane

CBSR can be constructed simply by replacing the Borel mass $M^{2}$ in Eqs. (1) and (2) by its complex generalization $\mathcal{M}^{2}$. As a further task, we have to choose the domain of complex Borel masses used in the MEM analyses. Because spectral functions obtained from MEM generally depend on the choice of the domain, the domain used for each channel in our analyses should be kept constant throughout the analyses at various temperatures, such that the temperature dependence of the spectral function can be extracted without artificial effects. In our analyses, the domain is fixed to the one used in vacuum for each channel. By considering the convergence of the OPE, the lower boundaries are determined as follows:

[TABLE]

where $d^{4}(\mathcal{M}^{2};T=0)$ is the whole dimension 4 term of the OPE in vacuum 222 This method of setting the lower boundaries was proposed in the first analysis using CBSR with MEM Araki et al. (2014)..

For the determination of the upper boundary, we employ an improved scheme as follows. In the originally proposed scheme Araki et al. (2014), we employed a simple circular boundary on the complex Borel mass plane, whose radius $M_{r}^{2}$ is treated as a free parameter as follows: $M^{2}<M^{2}_{r}$, where $\mathcal{M}^{2}=M^{2}\mathrm{e}^{i\theta}$. However, it includes the region $\theta\sim\frac{\pi}{2}$, where the damping by the kernels becomes weaker, as shown in the definitions

[TABLE]

[TABLE]

Because, in such a region, the integrals over the spectral function include large contributions from the continuum, it is natural to set a lower limit for the power of the exponential in these kernels as follows:

[TABLE]

Such a limit can adequately control the contribution from higher energy regions of the spectral function. It leads to a curve boundary on the complex Borel mass plane, as shown by the dashed line in Fig. 1. We choose the critical value $\tilde{r}_{c}$ such that the statistical significances of both the first and second peaks in vacuum become best. As a result, $\tilde{r}_{c}$ is equal to $1/2.2\,\mathrm{GeV}^{-2}$ and $1/2.0\,\mathrm{GeV}^{-2}$ for the V and PS channels, respectively. The actual domains for the V and PS channels determined by the above conditions, Eq. (8) and Eq. (11), are shown in Fig. 1 as shaded regions, in which the discretized complex Borel masses are used as input for our MEM analyses.

The default model, which is another MEM input (see e.g. Ref. Asakawa et al. (2001)), is chosen as a constant fixed to the asymptotic value of the spectral function computed by perturbation theory Morita and Lee (2010): the value of $\frac{1}{\pi}\mathrm{Im}\Pi_{pert.}^{V,PS}(\omega^{2})$ at $\omega=\infty$ for the V channel and $\omega=10\,\mathrm{GeV}$ for the PS channel 333Because $\mathrm{Im}\Pi_{pert.}^{PS}(\infty)$ diverges, we are forced to choose a finite energy to determine the “asymptotic value” of the spectral function. We checked that the obtained peaks do not strongly depend on the choice of this parameter.. The values of the other input parameters used in our analyses are summarized in Table 1.

III Analysis results

III.1 In vacuum

The spectral functions of both V and PS channels in vacuum obtained by our analyses are shown in Fig. 2. The estimated errors of the MEM results are shown by the three horizontal lines at each peak. It is observed that in both channels two clear peaks are generated. Both of them are statistically significant because their error bars lie between top and bottom of the peaks Asakawa et al. (2001); Gubler and Oka (2010). Note that, although a small third peak is also seen in both channels, it could be an artificial peak because it is not statistically significant. The positions of the first and second peaks agree with the experimental values with a precision of the order $50-150\,\mathrm{MeV}$ as shown in Table 2. Thus they are considered to correspond to the physical states, $J/\psi$, $\psi^{\prime}$, $\eta_{c}$ and $\eta_{c}^{\prime}$, respectively. It is interesting to note that all our obtained masses are lower compared to the experimental values. It is possible that this situation could improve once higher order $\alpha_{s}$ corrections to the Wilson coefficients are included in the analysis. As such terms have no influence on the temperature dependence of the gluon condensates, it is however not expected that they will qualitatively change our results about the temperature dependence of the peaks.

Let us note that, in an earlier study using conventional Borel sum rules with MEM Gubler et al. (2011), only one statistically significant peak was extracted for both channels, and it is widely distributed in the energy region between about $2.9\,\mathrm{GeV}$ and $3.6\,\mathrm{GeV}$. This wide peak likely corresponds to a combination of the first and second peaks obtained in our analysis. In other words, by the help of the higher resolution of CBSR, we are able to separate the degenerated of the previous work into two distinct peaks.

III.2 At finite temperature

Let us see how the peaks are deformed by finite temperature effects. Fig. 3 shows the spectral functions in the temperature range between $0.9\,T_{c}$ and $1.2\,T_{c}$ and also in vacuum for comparison. We observe in the figure that all the peaks tend to disappear with increasing temperature, that is, the charmonia melt by the temperature effects. In both channels, the spectral functions seem to be almost unchanged in the temperature range between $1.1\,T_{c}$ and $1.2\,T_{c}$. This observation may be interpreted as the fact that the complete melting of both states up to $1.1\,T_{c}$. On the other hand, the drastic changes of the spectral functions occur around at $1.0\,T_{c}$ for both peaks. This is understood to be caused by a sudden change of the $G_{0}(T)$ and $G_{2}(T)$ at about $1.0T_{c}$ (see e.g. Ref. Morita and Lee (2008a)).

Since we are most interested in comparing the melting behaviors of ground and excited states, we next examine them more carefully in the temperature region around $1.0\,T_{c}$. The bottom panels of Fig. 3 hence show spectral functions at temperatures between $0.96\,T_{c}$ and $1.04\,T_{c}$ with a interval of $0.02\,T_{c}$, while the energy range of the plot is limited to $\omega=2-5\,\mathrm{GeV}$. In this figure, however, we do not see any considerable difference between the ground and excited peaks: Both peaks seem to melt almost similarly with increasing temperature. On the other hand, in view of the MEM error estimation, we find that the second peak loses its statistical significance at a slightly lower temperature than that of the first peak, as shown in Table 3.

We note that the temperatures in Table 3 do not necessarily correspond to the original melting temperatures, but only indicate their “lower limits”. In Fig. 4, we show typical behaviors of error bars at finite temperature. At lower temperature $T<T_{\mathrm{low}}$, the peak has statistical significance, so that we can clearly conclude that it survives. At $T=T_{\mathrm{low}}$, the lower error bar of the peak reaches “the dip” at the high-energy side of the peak. At higher temperature $T>T_{\mathrm{low}}$, the statistical significance of the peak is lost, so that we cannot conclude anything about its existence. Thus, all we can say is that the physical peak survives at least up to $T=T_{\mathrm{low}}$, and $T_{\mathrm{low}}$ merely indicates a lower limit of the original melting temperature. Therefore, the original melting temperatures might become larger than the lower limits in Table 3.

As shown in Fig. 3 it is worthwhile to note that the positions of both peaks shift to the lower energy side of the order of $50\,\mathrm{MeV}$ before they lose their statistical significances. Such a shift was already obtained in earlier studies using Borel sum rules with/without MEM only for the ground states Morita and Lee (2008a, b); Song et al. (2009); Dominguez et al. (2010); Morita and Lee (2010); Dominguez et al. (2011, 2013); Kim et al. (2016); Kim and Lee (2017); Gubler et al. (2011); Suzuki et al. (2013) 444We note that the scattering term used in Ref. Dominguez et al. (2010), as shown in Appendix A of the reference, has an error.. In our analyses by CBSR with MEM, the same behavior is observed also for the excited states.

In MEM analyses, the shape of the obtained spectral function generally depends on the error of the input parameters. It is known that if the error grows larger, the spectral function tends to approach the default model, resulting in an unphysical suppression of a potential peak. In our QCD sum rule analysis with MEM, the error of the whole OPE becomes larger with increasing temperature through the growing uncertainties of $G_{0}(T)$ and $G_{2}(T)$. To test how large such unphysical effects are, we repeated the same analyses at finite temperature by keeping the error fixed to $T=0$ 555A similar analysis was performed also in Ref. Gubler et al. (2011).. By comparing such a test analysis and the original one at each $T$, we confirmed that the behavior of the respective spectral functions is not changed. Therefore, we can conclude that the deformation of spectral function obtained by our analyses is caused not by unphysical effects from the errors but by physically meaningful ones caused by the changing condensate values.

As a more quantitative analysis, we could fit the obtained spectral functions by a functional form with some fitting parameters (e.g. two peaks and continuum), as performed in Ref. Suzuki et al. (2013). However we did not obtain a stable solution for the temperature dependence of residues due to the appearance of local minima around $T_{c}$, where the spectral structures are more complex than the form of one peak and continuum as the previous study Suzuki et al. (2013). Thus we do not discuss residues at finite temperature in this paper.

IV Conclusion and outlook

We have applied CBSR with MEM to the vector and pseudoscalar channels of charmonium spectra in vacuum and at finite temperature. With the help of the higher resolution of CBSR, two statistically significant peaks in both channels are extracted, which were not resolved in the earlier Borel sum rules with MEM analysis of Ref. Gubler et al. (2011). Their peak positions agree well with the experimental values of $J/\psi$, $\psi^{\prime}$, $\eta_{c}$ and $\eta_{c}^{\prime}$ within $50-150\,\mathrm{MeV}$. By introducing finite temperature effects through two dimension-4 gluon condensates, we have observed that all the obtained peaks are deformed to gradually disappear as the temperature increases. We have confirmed that this is not an artificial MEM effect induced by the increasing error. They completely melt at the temperature of $T=1.1\,T_{c}$ since the spectral function does not largely change between $1.1\,T_{c}$ and $1.2\,T_{c}$ compared with drastic change around at $1.0\,T_{c}$. In both the channels, the first and second peaks seem to melt almost simultaneously. It is interesting to compare our predictions for excited states with recent experimental Khachatryan et al. (2014); Sirunyan et al. (2017); Adam et al. (2016a, b); Adare et al. (2017); Aaij et al. (2016) or theoretical results Aarts et al. (2011); Burnier et al. (2015); Braga et al. (2016).

To further improve the QCD sum rules, we can include higher dimensional terms in the OPE and their temperature dependences. For example, recently, the temperature dependences of the dimension-6 gluon condensates were phenomenologically estimated in Refs. Kim et al. (2016); Kim and Lee (2017). In the future, it would be worthwhile to precisely determine such higher dimensional condensates from lattice QCD simulations and to analyze QCD sum rules including these effects.

Finally, we comment on the possibility of studying bottomonium channels within the same approach. Recent experimental results from LHC show that the excited states, $\Upsilon(2S)$ and $\Upsilon(3S)$, are more strongly suppressed than the ground state $\Upsilon(1S)$ Chatrchyan et al. (2011, 2012, 2014); Khachatryan et al. (2017). Although we have tried to reproduce the bottomonium spectra using CBSR with MEM, we just obtained one single peak, generated from the overlap of the three original peaks, meaning that the resolution of MEM was not good enough to disentangle the three states for the bottomonium case. This is the same situation as usual Borel sum rules with MEM in our previous study Suzuki et al. (2013). The reason for this difficulty is related to the fact that the ratio of the mass differences between $\Upsilon(1S)$, $(2S)$, and $(3S)$ to the typical bottomonium mass scale $\sim 10\mathrm{GeV}$ is too small to separate these states by our approach, which is different from the analyses of charmonia with the mass difference of $m_{2S}-m_{1S}>500$ MeV and the energy scale $\sim 3\mathrm{GeV}$. The task of further improving the resolution of CBSR and finally separating the bottomonium states is left for future studies.

Acknowledgements.

We are grateful to Kenji Morita for providing us the gluon condensate data extracted from lattice QCD simulations. K.J.A. was supported by Grant-in-Aid for JSPS Fellows from Japan Society for the Promotion of Science (JSPS) (No. 15J11897). The research of P.G. is supported by the Mext-Supported Program for the Strategic Foundation at Private Universities, “Topological Science” (No. S1511006). This work is supported by Grants-in-Aid for Scientific Research from JSPS [Grant No. JP25247036(A)].