Quarkonia at T > 0 and lattice QCD

We report here on recent progress in the determination of S-wave and P-wave heavy-quarkonium states at finite temperature. Our results are based on the combination of effective field theories with numerical lattice QCD simulations. These non-perturbative tools allow us to compute the heavy-quarkonium in-medium spectral functions, from which we in turn determine the melting temperatures of individual states and estimate phenomenologically relevant observables, such as the ψ′ to J/ψ ratio in heavy-ion collisions.


Introduction
The bound states of a heavy quark and anti-quark (cc, bb), christened heavy-quarkonium (charmonium, bottomonium), have matured into a high precision probe in the study of heavy-ion collisions. Their constituents are produced early on in the partonic stages of a collision and subsequently sample the full evolution of the quark-gluon plasma (QGP) created in the collision center. Both the experimental programs at RHIC and LHC have over the last decade accumulated an unprecedented amount of high-precision data on these states. Among the highlights from run1 at the LHC [1] are the observation of relative excited states suppression in the bottomonium S-wave channel by the CMS collaboration, as well as the finding by the ALICE collaboration that the yields for the charmonium ground state J/ψ are replenished at √ s NN = 2.76TeV, compared to those measured at RHIC. The availability of such detailed measurements now urges theory more than ever to provide first-principle insight into the physics of heavy QQ in heavy-ion collisions. When we wish to start from first principles, naturally some idealization of the actual situation encountered in experiment is necessary. Here our main assumption is that of full kinetic equilibration of the heavy quarks with their surrounding, which is furthermore taken to be a thermal medium of quarks and gluons at a fixed temperature T . Since the temperatures encountered in current heavy-ion collision are not higher than 3−4×T PC with T PC = 155±5MeV [2,3] the chiral crossover temperature, we will use non-perturbative lattice QCD simulations to describe the medium degrees of freedom.
In such a static setup, we may compute the thermal spectral functions of heavy quarkonium, which provide us with vital information about their in-medium behavior. At T = 0 bound states appear as almost delta-like peaks, positioned at the mass of the state and well separated from the open heavyflavor threshold, at which the spectral function raises to a continuum. Their binding energy, both in vacuum and in-medium, is defined from the distance between the peak position and this continuum. At finite temperature these peaks may shift and thermally broaden and will eventually dissolve into a e-mail: rothkopf@thphys.uni-heidelberg.de the continuum. Besides the position and width of peaked structures, their area plays an important role too, as it is linked with the decay of the corresponding in-medium state either into dileptons or light hadrons [4]. I.e. using the position, width and area of the in-medium spectra, we will set out to estimate what consequences in-medium modification will have on measured yields in heavy-ion collisions.

Quarkonium in-medium spectral functions
For theory, heavy-quarkonium represents a unique probe, due to the inherent separation of scales between the heavy-quark rest mass m Q and e.g. the temperature T/m Q 1 or the characteristic scale of QCD Λ QCD /m Q 1. In the presence of such a separation, we may use the small ratios as expansion parameters to systematically go over from the language of quantum chromodynamics (QCD) to a non-relativistic description of the heavy quarks in terms of effective field theories (EFT) [5]. We will hence compute the in-medium quarkonium spectral functions in the following using two complementary EFT approaches combined with lattice QCD simulations: The first and direct one is based on non-relativistic QCD (NRQCD) [6,7], which treats realistic heavy-quarkonium with a finite mass. However with the currently available simulation data its resolution is still limited. This line of study is also pursued by the FASTSUM collaboration in Refs. [8][9][10][11]. The second approach is an indirect one, based on the EFT potential NRQCD (pNRQCD) (for a perturbative study see [12,13]), where one first computes the in-medium potential between the QQ and in turn solves a Schrödinger equation for the in-medium spectral function. And while this approach is not limited in resolution, the underlying T > 0 potential we use at the moment does not yet include finite velocity or spin dependent corrections.
Both approaches face the same challenge that at some stage they need to extract dynamical information in the form of spectral functions from lattice QCD simulations, which are carried out in an unphysical Euclidean time. The simulated correlation functions D(τ i ) = D i are related to the spectral functions via a Laplace-type integral transform, which needs to be inverted Since we need to discretize the spectrum at many more frequencies, than we have simulation data available and each D i comes with a finite error, the inversion is inherently ill-defined. I.e. a χ 2fit of the ρ l 's would lead to an infinite number of degenerate solutions that all reproduce D i within uncertainties. One possibility to nevertheless give meaning to the inversion is to use Bayes theorem [14], which provides a systematic prescription how to regularize the naive χ 2 fit by incorporating additional knowledge (I) on the spectra, e.g. positive definiteness or smoothness conditions. This prior knowledge is usually encoded in a function m l called the default model, which represents the extremum of the regulator functional, the prior probability P Once the regulator and m is specified, one ends up with a numerical optimization problem to determine the most probable spectrum given simulation data and prior information. In the following we will use two different implementation on the market, the well known Maximum Entropy Method (MEM) [15] and a more recent Bayesian Reconstruction (BR) method [16]. They differ both in the form of the regulator and the implementation of the numerical search for ρ Bayes . Over the last two years we have  gained a much better understanding of the very different systematic uncertainties of the two methods: the MEM has been found to be prone to over-smoothing, in particular if only a relatively small number of datapoints is available, while the BR method can introduce ringing artifacts that may mimic peak features not actually present in the simulation data. Note that as long as we are not in the "Bayesian continuum limit", i.e. N d → ∞ and ∆D/D → 0 two different implementation of the Bayesian strategy may provide different answers and only in that limit will agree. Unfortunately due to the strong loss of information due to the integral Kernel in (1) approaching this limit is exponentially hard. These challenges will be apparent in the direct analysis of spectral functions using NRQCD on the lattice in the following subsection.

Direct determination based on lattice NRQCD
Instead of simulating the light and heavy d.o.f. on the same space-time lattice, which would require too fine of a lattice spacing, we here [17] use the separation of scales T/m Q 1, λ QCD /m Q 1 to separate the physics of the thermal medium from that of the heavy quarks in the following way. In this lattice variant of the EFT NRQCD, the light quarks and gluons are treated with a standard lattice QCD simulation at T > 0, while the heavy quarks are described by non-relativistic Paulispinors propagating in the background of the thermal fields. At this point no modelling is involved, as NRQCD is systematically derived from the QCD Lagrangian by an expansion in 1/m Q a with a being the lattice spacing. Current implementations include up to fourth order corrections, corresponding to O(v 4 ) in the language of continuum NRQCD. For the medium we utilize state-of-the-art lattices by the HotQCD collaboration [18] on 48 3 × 12 grids featuring dynamical u, d and s quarks based on the HISQ formulation with almost physical pion mass m π = 161MeV. These simulations span temperatures T ∈ [140 . . One then computes the propagation of a single heavy quark in the medium background, and combines two of these propagators into a heavy quarkonium correlator. Projected to a definite quantum number by the insertion of appropriate vertex operators one arrives at a heavy quarkonium correlation function in Euclidean time. This quantity can already tell us about overall in-medium modification if we divide its value at T > 0 by that in vacuum, as shown in Fig. 1 for Upsilon (left) and χ b (right).
If the ratio is unity within statistical errors, no in-medium modification is present. Around T C we find only very mild deviations from one, albeit there is a hint for a non-monotonic behavior present.   At higher temperatures a characteristic upward bending occurs but even at the highest T = 407MeV the maximum departure from unity is 1.75% for Upsilon and 6.5% for χ b . The differences between the two cases however are a first sign that in-medium modification is actually hierarchically ordered with the vacuum binding energy of the ground state of a particular channel.
From the in-medium correlators we may now reconstruct the spectral functions using both the BR method and the MEM. We use a flat default model m in the following and discretize frequencies with N ω = 3000 leading us to the results for the Upsilon channel shown in Fig. 2. Note that due to the small number of simulation datapoints N d = 12 available only the ground state peak structure is captured robustly, all other higher lying features must be attributed to numerical artifacts.
We find that the BR method shows a well pronounced lowest lying structure at all temperatures up to T = 407MeV, i.e. the amplitude of the lowest peak is at least a factor 5 larger than the amplitude of the neighboring artifacts. The MEM on the other hand shows only washed out features above T ≈ 333MeV. Similar MEM results have been obtained by the FASTSUM collaboration in [10].
To interpret our results we need to consider the systematic uncertainties of the two different methods. The possibility for over-smoothing in the MEM and ringing in the BR method. I.e. one is bracketing the actual temperature for the disappearance of a state from above and below with the two methods. In the Bayesian continuum limit these temperatures would become the same. Now let us continue to the P-wave states given in Fig. 3. The BR method results shown on the left, at first sight, seem to indicate that again up to T = 407MeV a well defined lowest lying peak is present, while the MEM, consistent with previous FASTSUM computations shows only washed out features above T ≈ 210MeV.
A simple inspection by eye is however not enough in case of the BR method. Since the amplitude of the lowest peak at T = 407MeV is already lower than that of the neighboring numerical artifact at higher frequencies, we need a more systematic way to determine whether that feature really represents a sign of a remnant bound state. Our strategy to do so is to compare with reconstructed non-interacting spectral functions. The actual free spectrum is known to be a continuous functions devoid of peaks. However if such spectrum is reconstructed from a small number of correlator points, akin to Gibbs ringing in the inverse Fourier transform we expect the BR method to also introduce wiggly behavior.  to the vacuum ground state at all temperatures, while (right) the MEM shows only washed out features above T ≈ 210MeV. As described in the text, we can identify that the apparent ground state feature at T = 407MeV in the BR method is a numerical artifact and thus both methods see the disappearance of χ b below T = 407MeV. Due to the different systematics of the two methods they however bracket the actual disappearance of the bound state peak from above (BR) and below (MEM). And indeed as shown in Fig. 4 the gray curve corresponding to the reconstructed free spectrum shows sizable ringing artifacts. In particular their amplitude is at least as strong as that of the apparent bound state signal in the interacting reconstructed spectrum. We conclude that the BR method does not give any indication of a bound state remnant at T = 407MeV and note again that by using the two different methods, BR and MEM, we are thus bracketing the actual disappearance of the χ b bound state from above and below. While this direct approach to in-medium quarkonium based on lattice NRQCD is promising the available simulation data needs to be improved to access also the physics of e.g. excited states in this framework. In the spirit of reaching the Bayesian continuum limit, we are thus currently working on increasing the statistics, while e.g. the FASTSUM collaboration is generating ensembles with a larger number of points N d along Euclidean time direction.

Indirect determination based on pNRQCD
A complementary but indirect approach to in-medium heavy quarkonium spectral functions is based on the effective field theory pNRQCD. In it the QQ two-body system is described in terms of nonrelativistic color singlet and color-octet wave-functions. The correlation functions of such wavefunctions evolve under a Schrödinger like equation of motion featuring a potential V QCD (r) that encodes the interactions between the constituents and the environment. As the EFT can be systematically derived from QCD this potential consists of a lowest order static contribution and may be amended by velocity and spin dependent corrections. Depending on the degree of scale separation also nonpotential effects may play a role but are not further considered here. The values of V QCD are obtained from the underlying microscopic field theory by the process of matching, where two correlation functions of equal physical content, one in pNRQCD, one in QCD are set equal at a certain scale. It has been shown that the in-medium potential can be related to a real-time quantity in thermal QCD, the rectangular Wilson loop Evaluating this definition in resummed perturbation theory at high temperature [19,20] Based on a general argument using the symmetries of the real-time Wilson loop we have shown [23] that if a well defined lowest lying spectral feature exists in ρ a potential can be meaningfully defined. The values of Re[V QCD ] and Im[V QCD ] are furthermore related to the position and width of that lowest lying spectral feature. And since the Wilson loop spectral function is actually much simpler in structure than the quarkonium spectral function we are able to reconstruct it already with currently available simulation data (Note that in order to avoid cusp divergencies, in practice we compute the Wilson line correlator in Coulomb gauge instead of Wilson loops). Here we use [24] lattice QCD simulations by the HotQCD collaboration with dynamical u, d and s quarks based on the asqtad formulation with heavier than physical pion mass m π ≈ 300MeV. The corresponding outcome is shown as filled points in Fig. 5 for Re[V QCD ] on the left and Im[V QCD ] on the right. The real part is shifted by hand in y for better readability and shows a clear, as well as smooth transition from a confining behavior close to T = 0 to a screened form above T = T C . In the left panel we show as light colored points the values obtained for the imaginary part (shifted in x), which is finite above T C . The collection of discrete points obtained so far however is not enough to compute the heavy quarkonium spectral function, to this end an analytic parametrization is required.
Recently we derived such an analytic parametrization [25] starting from the realization that at T ≈ 0 the values of Re[V QCD ] we obtain on the lattice are very well described by a simple Cornell ansatz V QCD T=0 (r) = V Cornell (r) = −α S /r + σr + c. α S denotes the strong coupling, σ the string tension and c is an arbitrary scale dependent shift. These three parameters in the following remain temperature independent and are determined via a fit to the T ≈ 0 datapoints on the left of Fig. 5. The corresponding best fit is given as the solid dark violet and dark blue line. Similar to the Gauss law for the Coulombic part of the Cornell potential one can establish a generalized Gauss law to also describe the linearly rising part.
In-medium effects are introduced in this setup by using the in-medium permittivity of a weakly coupled gas of quarks and gluons, taken from hard-thermal loop resummed perturbation theory. The idea is that the non-perturbative effects of inter-quark binding and confinement are encoded in the Cornell form of the T = 0 potential and that the HTL contributions to the medium will as a first step   Fig. 6, and are consistent with zero below T C while quickly switching to finite values above. The red solid line shows a fit, which approaches the perturbative HTL value of m D at very high T T C . To compute realistic heavy quarkonium spectra of course we would have to use continuum extrapolated values for the potential. As these are not yet available but still work in progress, we will use the following strategy. At T = 0 we fit the values of the Cornell potential parameters so as to reproduce the spectrum of bottomonium bound states listed in the PDG. The in-medium modification is then imprinted via the Gauss-law parametrization and the continuum corrected values of m D obtained on the lattice.  With all these ingredients set, we can now compute the in-medium spectral functions by solving the Schrödinger equation for the so called forward correlator D > (t, r), which is in essence the unequal time correlator of singlet wavefunctions [13]. Note that this is not the same as solving the Schrödinger equation for an individual realization of the wavefunction. The sought after quarkonium spectrum is related to this correlator via a Fourier transform and an appropriate limiting procedure lim r→0 dte −iωt D(t, r) = ρ(ω). The results for the S-wave channel [26] are shown in Fig. 7 and for the P-wave [27] in Fig. 8, the left panel contains the spectra for bottomonium, the right for charmonium respectively.
We observe characteristic changes in all the channels, i.e. a broadening and shifting of individual peaks to lower masses as temperature is increased. The strength of the modification is clearly ordered according to the vacuum binding energy of the individual states. Besides the bound state features also the position of the heavy-flavor threshold moves to lower energies, since the real-part of the potential becomes more and more weakened as T rises. Eventually this continuum engulfs the bound state remnants and they disappear. This however is a gradual process and one needs to set a quantitative criterion at which T we declare a state melted. With the information from the spectra available we may use the popular choice of Γ(T ) = E bind (T ) (see e.g. [19]), i.e. the point at which the in-medium width equals its in-medium binding energy. In a naive wavefunction picture it correspond a dampening by a factor 1/e after one oscillation. The corresponding melting temperatures we obtain for the S-wave and P-wave are given below  (2) Melting temperatures are interesting in their own right, however they are not directly measurable in experiment. Therefore in the next section we wish to discuss our proposal how to proceed towards extracting phenomenologically relevant observables, such as the ψ to J/ψ ratio from the in-medium spectra.

In-medium quarkonium phenomenology from spectral functions
To connect to experimental results, we need to keep two caveats in mind. First, one does not measure in-medium dilepton emission of heavy quarkonium in experiment but instead the decay of vacuum states long after the QGP has ceased to exist. Therefore comparing the in-medium spectra directly to measured spectra is not admissible. On the other hand our assumption of full kinetic equilibration is, if at all, only applicable to charmonium at low p T and mid-rapidity. Therefore in the following we concentrate the lighter of the two flavors.
Then we may ask what are meaningful observables for charmonium? Take the nuclear modification factor R AA of J/ψ for example. Several models based on different physics assumptions do reproduce this quantity fairly well, making it difficult to pinpoint the relevant processes unambiguously. Therefore it has been proposed to investigate more discriminatory observables, such as the ψ to J/ψ ratio, for which prediction by transport models or the statistical model of hadronization differ significantly. Here we wish to estimate this ratio using the in-medium spectra of fully equilibrated charmonium. T>0 spectra Figure 9. Our estimate for the ψ /J/ψ ratio (orange) as well as that from the statistical model of hadronization (pink) and a recent first measurements by ALICE (adapted from [1]).
The additional assumption we have to make here is that of an instantaneous freezeout, i.e. similar to the statistical model of hadronization, at T C = 155MeV we let all in-medium state go over their vacuum counterparts. We do so in the following way: the question we ask is how many vacuum states does the in-medium spectral peak in the charmonium spectral function correspond to? The answer we give in units of dilepton emission, which is related to the weighted area R ll ∝ d p 0 d 3 p ρ(P) P 2 n B (p 0 ) under each spectral peak. I.e. we compute R ll for both the inmedium and the T = 0 spectrum and divide the two, leading us to an estimate of the number of states produced. If we do so for both J/ψ and ψ and divide the obtained values the outcome is Compared to a recent first measurement of the ratio by the ALICE collaboration shown in Fig. 9, we find that our result is slightly larger than the prediction by the statistical model of hadronization (SHM) but still well within the limits given by ALICE. While the difference to the SHM is larger than the estimated uncertainties it is not unfathomable that two approaches that are based on a full kinetic equilibration of the heavy quarks do show such a similar outcome.

Conclusion
The study of in-medium properties of heavy quarkonium from first principles has progressed significantly over the last years. The maturation of effective field theories, the arrival of realistic lattice simulations of QCD at finite temperature with almost physical pion mass and the development of novel Bayesian approaches to spectral function reconstruction have all contributed to further the computation of quarkonium in-medium spectral functions. Here we reported on two complementary approaches combining lattice QCD with the effective field theories NRQCD and pNRQCD. The former, direct approach, while describing realistic heavy-quarkonium with finite mass is still limited in resolution due to the small number of simulation datapoints available. The measured correlators indicated that overall in-medium modification proceeds hierarchically ordered according to the vacuum binding energy of the ground state in each channel. Using two different Bayesian spectral reconstructions (BR) and (MEM) we reconstructed the corresponding in-medium spectra and bracketed the melting temperatures of bottomonium S-wave and P-wave ground states from above and below. The systematic uncertainties of the different reconstruction methods were discussed. The latter, indirect approach, as a first step required the computation of the in-medium heavy-quark potential, for which we showed our most recent results for full lattice QCD. Based on a generalized Gauss law in which the in-medium modification of Re[V QCD ] and Im[V QCD ] is described by a single temperature dependent parameter m D , we managed to describe our lattice data for the potential. Using the continuum corrected in-medium potential we then solved a Schrödinger equation to obtain the inmedium quarkonium spectra. They showed a hierarchical modification of the former vacuum states according to their binding energy. The peaks broaden and shift to lower masses, before vanishing into the continuum, which at the same time also moves to lower and lower energies. Using the criterion of the in-medium width being equal to the in-medium binding energy we defined and evaluated the melting temperatures of the individual states.
With the additional assumption of an instantaneous freezeout we estimated the ψ to J/ψ number ratio obtaining a value slightly larger than the estimate from the statistical model of hadronization.
The author thanks Y. Burnier, O. Kaczmarek, S. Kim and P. Petreczky for the fruitful collaboration in the presented studies. Calculations for [26] were performed on the in-house cluster at the ITP in Heidelberg and the SuperB cluster at EPFL. This work is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)".