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

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.


I. INTRODUCTION
The particle consisting of a heavy quark and an antiquark QQ, "quarkonium", has been a suitable target to study dynamics of QCD at short distance due to its large mass [1]. 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 [2]. 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 QQ potential [2][3][4][5][6][7][8][9][10][11][12] or the disappearance of peaks in the spectral function by the finite temperature effects. The temperatures at which the peaks completely disappear, "melting temperatures" 1 , 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) [13] 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/ψ and η c ) survive above the temperature of 1.5 T c [14] (cf. [15][16][17][18][19][20][21]). Similarly, QCD sum rules [22,23] can study quarkonium suppression by incorporating finite temperature effects through the QCD * k.araki@th.phys.titech.ac.jp † k.suzuki.2010@th.phys.titech.ac.jp ‡ pgubler@riken.jp § oka@th.phys.titech.ac.jp 1 Actually, 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.
In the recent heavy-ion collision experiments, the suppression of charmonium excited states (ψ or ψ(2S)) was also observed [35][36][37][38][39][40]. 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 [2][3][4][5][6][7][8][9][10]. 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 [34].
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 [33,34,[41][42][43][44][45][46], 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 [47]. In Ref. [47], the excited φ 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. [33,34]. For specifying the domain of the complex Borel plane to be used in the MEM analysis, we propose an updated criterion.

A. 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Π V (q 2 ) ≡ Π V,µ µ (q)/(−3q 2 ) andΠ P ≡ Π 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: where ν is a dimensionless parameter defined as ν = 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: Let us explain each part of Eq. (1). α s (ν) is a strong running coupling constant evaluated at the Borel mass scale M . The functions, A J , a J , b J and c J are the Wilson coefficients. Their explicit form can be found in Ref. [28]. The first and second terms come from the usual perturbative result up to the first order in α s . The φ(T ) s are the condensates defined as 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: where the expression O T means the Boltzmann average defined as O T = Tr(e −H/T O)/Tr(e −H/T ). u is the four velocity vector of the medium whose norm is unity.

B. 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 ∼ 1, GeV. To proceed further, we therefore need to calculate the temperature dependences of the gluon condensates. In our formulation, the strategy proposed in Refs. [24,25] is employed, in which we express the two gluon condensates as follows: where , p and α 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 [48,49], 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 MeV, while the full QCD result leads to the cross-over temperature being 145 − 165 MeV [50,51]. 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 qq and qγ µ D ν q , which will hence modify Eqs. (6) and (7). Once the temperature dependences of all these operators (and α s ) are known in full QCD, it will become possible to perform a QCD sum rule calculation that goes beyond the quenched approximation.
C. 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 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: where d 4 (M 2 ; T = 0) is the whole dimension 4 term of the OPE in vacuum 2 . For the determination of the upper boundary, we employ an improved scheme as follows. In the originally proposed scheme [47], we employed a simple circular boundary on the complex Borel mass plane, whose radius M 2 r is treated as a free parameter as follows: M 2 < M 2 r , where M 2 = M 2 e iθ . However, it includes the region θ ∼ π 2 , where the damping by the kernels becomes weaker, as 2 This method of setting the lower boundaries was proposed in the first analysis using CBSR with MEM [47].m c(mc) 1.273 ± 0.006 [52] αs π G 2 0.012 ± 0.0036 GeV 4 [23, 53] ΛQCD 0.213 ± 0.008 [54]  shown in the definitions 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: cos θ M 2 >r c .
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 valuer c such that the statistical significances of both the first and second peaks in vacuum become best. As a result,r c is equal to 1/2.2 GeV −2 and 1/2.0 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. [13]), is chosen as a constant fixed to the asymptotic value of the spectral function computed by perturbation theory [28]: the value of 1 π ImΠ V,P S pert. (ω 2 ) at ω = ∞ for the V channel and ω = 10 GeV for the PS channel 3 . The values of the other input parameters used in our analyses are summarized in Table I. 3 Because ImΠ P S pert. (∞) 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.

A. 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 [13,41]. 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 MeV as shown in Table II. Thus they are considered to correspond to the physical states, J/ψ, ψ , η c and η c , 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 α 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 [33], only one statistically significant peak was extracted for both channels, and it is widely distributed in the energy region between about 2.9 GeV and 3.6 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.

B. 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. [24]).
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 ω = 2 − 5 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 III.
We note that the temperatures in Table III 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 low , the peak has statistical significance, so that we can clearly conclude that it survives. At T = T low , the lower error bar of the peak reaches "the dip" at the high-energy side of the peak. At higher temperature T > T 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 low , and T 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 III. 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 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 [24][25][26][27][28][29][30][31][32][33][34] 4 . 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  Table III. Right (T > T low ): Statistically nonsignificant peak.
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 5 . 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. [34]. 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 [34]. 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. [33]. Their peak positions agree well with the experimental values of J/ψ, ψ , η c and η c within 50−150 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 [35][36][37][38][39][40] or theoretical results [55][56][57].
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. [31,32]. 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, Υ(2S) and Υ(3S), are more strongly suppressed than the ground state Υ(1S) [58][59][60][61]. 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 [34]. The reason for this difficulty is related to the fact that the ratio of the mass differences between Υ(1S), (2S), and (3S) to the typical bottomonium mass scale ∼ 10GeV 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 ∼ 3GeV. The task of further improving the resolution of CBSR and finally separating the bottomonium states is left for future studies.