Quarkonium at finite temperature: Towards realistic phenomenology from first principles

# Quarkonium at finite temperature: Towards realistic phenomenology from first principles

Yannis Burnier, b    Olaf Kaczmarek c    and Alexander Rothkopf Institute of Theoretical Physics, EPFL, CH-1015 Lausanne, SwitzerlandFakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, GermanyInstitute for Theoretical Physics, Heidelberg University, Philosophenweg 16, 69120 Heidelberg, Germany
August 18, 2019
###### Abstract

We present the finite temperature spectra of both bottomonium and charmonium, obtained from a consistent lattice QCD based potential picture. Starting point is the complex in-medium potential extracted on full QCD lattices with dynamical u,d and s quarks, generated by the HotQCD collaboration. Using the generalized Gauss law approach, vetted in a previous study on quenched QCD, we fit with a single temperature dependent parameter , the Debye screening mass, and confirm the up to now tentative values of . The obtained analytic expression for the complex potential allows us to compute quarkonium spectral functions by solving an appropriate Schrödinger equation. These spectra exhibit thermal widths, which are free from the resolution artifacts that plague direct reconstructions from Euclidean correlators using Bayesian methods. In the present adiabatic setting, we find clear evidence for sequential melting and derive melting temperatures for the different bound states. Quarkonium is gradually weakened by both screening () and scattering () effects that in combination lead to a shift of their in-medium spectral features to smaller frequencies, contrary to the mass gain of elementary particles at finite temperature.

## 1 Introduction

Heavy quarkonium, the bound states of a heavy quark and anti-quark ( or ), are a unique tool to simplify the complexity inherent in the physics of the strong interactions. Instead of having to deploy the full force of quantum field theory, we may consider a Schrödinger equation with an effective potential to describe the dynamics of these bound states in vacuum and in-medium. At zero temperature the two main characteristics of QCD, asymptotic freedom and confinement manifest themselves clearly in the form of a Coulombic behavior with a running coupling at small distances and a non-perturbative linear rise at large distances Koma:2006si . I.e. we can learn about key features of the strong interactions by inspecting this simple potential. Vice versa, from the knowledge of the potential, we can reproduce even quantitatively a significant part of vacuum quarkonium physics (e.g. the bound state spectrum below the open heavy quark threshold Quigg:1979vr ; Eichten:1995ch ).

The question of how such a potential arises from the microscopic theory of QCD has been answered in detail by advances made in the research on the effective field theories (EFT) NRQCD and pNRQCD. The latter has been introduced at first at in Brambilla:1999xf and investigated in the non-perturbative regime in Brambilla:2000gk ; Pineda:2000sz , with both cases being reviewed in detail in Brambilla:2004jw . The generalization to finite temperature in a perturbative setting followed in Brambilla:2008cx .

In the presence of a separation of scales, the combination of NRQCD and pNRQCD offers a systematic power counting prescription in the heavy quark velocity , to setup a simplified description of the bound state evolution at soft energy scales () in terms of singlet and color octet wavefunctions. I.e. we do not have to explicitly describe the physics at the much higher hard scale ()111The strength of the EFT approach lies in the fact that any residual influence of the physics of the hard scale can be systematically included via the matching of Wilson coefficients, which goes beyond the ability of direct methods, such as e.g. the historic Wilson loop approach Barchielli:1986zs . And indeed the heavy quark rest mass (in this work we use GeV and GeV) and the characteristic scale of quantum fluctuations in QCD, usually denoted by MeV are far enough apart to warrant a quantitatively reliable potential description in vacuum.

Effective field theory has furthermore reminded us that the potential between two infinitely heavy quarks is an inherently Minkowski-time quantity. Starting out from fundamental QCD, the evolution of a can be described by the thermally averaged unequal-time point-split meson-meson correlator

 D>(r,t)=⟨M(x,y,t)M†(x0,y0,0)⟩, (1)

with the meson operator defined as . refers to the Dirac matrices and denotes a straight Wilson line connecting and . The relative coordinate enters as .

In the limit of infinite quark mass, where the static constituents are truly test particles, (1) can be identified with the medium averaged unequal time correlation function of quarkonium singlet wavefunctions in the language of the EFT. It is this wavefunction correlator for which a Schroedinger equation and thus an in-medium potential will be established. For static quarks the spatial separation becomes an external parameter of the theory and their evolution traces out a rectangle over time. Formally it has been shown that eq.(1) reduces to the rectangular real-time Wilson loop Brambilla:2004jw

 W□(r,t)=⟨Tr(exp[−ig∫□dxμAaμTa])⟩. (2)

As such, this object still contains the physics of all scales: hard, soft and ultra-soft, i.e. in general it fulfills an equation of motion Laine:2007qy

 i∂tW□(r,t)=Φ(r,t)W□(r,t) (3)

with being a time- and space dependent complex function.

In case that a potential picture for the static system is applicable, the function has to asymptote to a time independent function . In turn it will dominate the evolution of at times much larger than the intrinsic scales of e.g. the gluons and light quarks in the medium. Formally we may then write

 V(r)=limt→∞i∂tW(t,r)W(t,r). (4)

This limit reflects the fact that in order to replace a retarded, i.e. gluon mediated interaction with an instantaneous potential, the gluons must have been exchanged between the heavy quarks multiple times.

In the presence of the additional scale we need to ascertain how reliable a description of realistic, i.e. finite-mass, quarkonium is in terms of the static potential alone. Two forms of so called relativistic corrections can adversely affect the accuracy of the lowest order approximation. The first kind remains fully within the potential picture, i.e. finite mass corrections to the static potential can become significant, as they are proportional to the relative velocity of the heavy quarkonium system. These corrections can be systematically computed (see e.g. Brambilla:2000gk ; Pineda:2000sz ; Koma:2006si ), in vacuum they are found to be small and we expect the same to be true at finite temperature. Their non-perturbative determination from finite temperature lattice QCD will be the aim of a future study.

The other contributions are so called non-potential effects, arising from physics that cannot be recast into the form of a time independent term in the Schrödinger equation. One way they can enter in the perturbative formulation of effective field theories if the ultra-soft gluons are still kept as explicit degrees of freedom (see e.g. Brambilla:2004jw ; Bazavov:2014soa ). In a non-perturbative setting, if non-potential effects become sizable, the late time evolution of the Wilson loop cannot be described by a simple time-independent potential . As will be discussed in the next section we find that at the temperatures and inter-quark distance investigated here, the size of such contributions remains insignificant.

Interestingly the definition of the potential from the Wilson loop in real-time coincides with the late limit in Euclidean time in the case of . This constitutes the basis for the successful extraction of the static zero temperature potential from lattice QCD simulations, which are carried out solely in an imaginary time setting. At finite temperature, the real-time definition of the potential does not change, however the imaginary time axis becomes compactified and its finite extend encodes vital physics information, i.e. the inverse temperature. Hence the straight forward connection between the late Minkowski time limit and the Euclidean Wilson loop at maximum is absent.

For more than two decades, the absence of an effective-field theory based definition for the the in-medium potential has led theorists to embrace model potentials that were defined directly from Euclidean time observables Nadkarni:1986cz , readily calculable in lattice QCD. In particular two quantities have gained popularity as model potentials, the color singlet free energies in Coulomb gauge

 eβF(1)(r)=⟨Tr[Ω(r)Ω†(0)]⟩CG,Ω(r)=exp[−ig∫β0dτAaμ(r,τ)Ta] (5)

defined from the correlator of Polyakov loops and a derived quantity the color singlet internal energies . While these quantities exhibit a behavior compatible with the expectations for e.g. Debye screening of the interaction between the heavy quarks in the deconfined phase, it could be shown that they do not match the potential for the quark anti-quark system at finite temperature Burnier:2009bk ; Brambilla:2010xn .

And indeed neither one of these fully real quantities by themselves can take the role of a static in-medium heavy quark potential, as we have learned from a ground breaking series of works starting with Laine et. al. Laine:2007qy ; Brambilla:2008cx in 2007. In their study the real-time definition (4) was evaluated in a resummed perturbative framework, called the hard-thermal loop approximation. This approach contains a gauge invariant resummation of an infinite number of Feynman diagrams and has been shown to capture many key features of QCD at high temperature reliably. As the Wilson loop in Minkowski time is an in general complex quantity, the authors observed that the potential too takes on complex values at finite temperature

 VHTL(r)=−~αs⎡⎣m\tiny\rm{D}+e−m\tiny\rm{D}rr+iTϕ(m\tiny\rm{D}r)⎤⎦+O(g4), (6)

with

 ϕ(x)=2∫∞0dzz(z2+1)2(1−sin(xz)xz). (7)

Here we absorb a factor in the definition of the coupling constant to connect to the conventions in the phenomenology literature. It was furthermore shown that the real-part of this complex potential itself does not coincide with the color-single free energies in hard-thermal loop resummed perturbation theory in Burnier:2009bk ; Brambilla:2010xn , even though the absolute deviations are comparatively small.

The fact that the in-medium potential is a complex quantity not only reflects a quantitative change but necessitates a qualitatively different perspective on the physics it describes. I.e. besides screening of the force between the heavy quarks (Debye screening) seen in the real-part, the effects of scattering of light medium degrees Laine:2007qy ; Brambilla:2013dpa with the heavy quarks (Landau damping) related to further weaken the bound state in the QCD medium. In fact, depending on the hierarchy of scales present, the imaginary part can be related to different phenomena, such as the breakup of a color singlet to an octet configuration Brambilla:2008cx and in turn to the dissociation of the system into gluons Brambilla:2011sg .

Only these effects taken together can give a consistent picture of heavy quark bound states at finite temperature. I.e. a simple estimate of the dissolution of a heavy quarkonium state solely on the basis of vanishing binding energy is not sufficient. This in turn has consequences for heavy quarkonium phenomenology in general, where e.g. the melting temperatures enter as input into transport model calculations.

We would like to stress that the in-medium potential discussed here does not directly govern the evolution of the actual wave-function of the heavy quarkonium system. By construction it instead enters in the Schrödinger equation of the real-time Wilson loop, which in the EFT language corresponds to the thermally averaged correlator of unequal time wavefunctions. The imaginary part of the potential hence describes the decay of this correlator over time, which does not directly imply the annihilation of the heavy-quarks. In fact, due to the non-relativistic approximation we operate under, and remain in the system forever. However, even if the norm of the quarkonium wavefunction is preserved in unitary time evolution, its correlation with the initial state can still decay with time, a phenomenon known as decoherence. It is an active area of research to answer the question how to connect the complex in-medium potential to the evolution of the bound state wave-function, which has not yet been answered in the effective field theory setting of pNRQCD. The concept of open-quantum systems (see e.g. Akamatsu:2011se ) has proven insightful in this regard.

One possible way to elucidate the physics encoded in the complex potential is to compute the spectral function of heavy quarkonium. This quantity is related to the heavy quark current-current correlator, which can be obtained from (1) by carefully taking the limit of vanishing point splitting. In section 3 we will hence solve an appropriate Schrödinger equation to compute . Once computed one can observe how the formerly delta-like bound state peaks present at broaden and shift as screening and scattering modifies the state with increasing temperature. Choosing a popular criterion of melting temperature Laine:2007qy ; Kakade:2015xua at the point where binding energy and spectral width coincide we can furthermore determine the point of dissolution for different bottomonium and charmonium states.

Changes in spectral structure are of particular interest as they are directly related to changes in the dilepton emission from quarkonium decay. The dilepton emission rate is given by a simple product of the in-medium spectral function with the Bose-Einstein factor McLerran:1984ay

 dRℓ¯ℓd4P=−Q2qα2e3π3P2nB(p0)ρ(P), (8)

where denotes the electric charge of the heavy quark considered in units of , the four momentum is and the finite mass of the leptons has been neglected (for a more detailed discussion of the above formula see ref. Laine:2013vma ). I.e. if a bound state or its remnant appears as well defined peaked feature in an in-medium spectral function, the area under such a structure informs us about the experimentally accessible dilepton emission of that state. Note that while the above relation is only applicable to the decay of a quarkonium state in a thermalized static plasma, it can nevertheless provide us with vital physics insight as we will see in sec. 4.1.

This paper is organized as follows. We start sec.  2 with a review on recent progress in the extraction of the values of the complex in-medium potential from Euclidean time lattice QCD simulations. In order to deploy the lattice potential in phenomenological applications, sec. 2.2 discusses the generalized Gauss law ansatz developed in ref. Burnier:2015nsa , which provides an analytic formula for and depending on a single temperature dependent parameter the Debye mass. Tuning we are able to reproduce the lattice values of and confirm the yet tentative values for . These values for are compared to perturbation theory. Section 3 is concerned with using the continuum corrected potential to investigate the phenomenology of in-medium heavy quarkonium. We will calculate the spectra for the bottomonium and charmonium S-wave channel and investigate their modification with increasing temperature. Besides giving the melting temperatures for different states, we will determine the to ratio at the chiral crossover and give a rough estimate for the suppression of bottomonium in a heavy-ion collision compared to . We close with a conclusion in sec.  5.

## 2 Lattice QCD potential and Debye screening mass

While the perturbative computation of the potential has contributed significantly to our understanding of the physics involved, it is not sufficient for the description of the experimentally relevant temperature regime around the phase transition. There the quark gluon plasma can indeed be considered strongly interacting, exemplified, e.g. by the large value of the trace anomaly Bazavov:2014pvz . This calls for an evaluation of the potential definition (4) in lattice QCD, which at first seems unfeasible, since direct access to real-time quantities, such as the Wilson loop (2) is not possible. Conceptual and technical advances in the extraction of real-time information from lattice QCD simulations over the last few years however have made such an evaluation possible, as we will discuss below.

### 2.1 The complex in-medium potential from lattice QCD

The real-time Wilson loop cannot be directly accessed in lattice QCD simulations, as they are performed in Euclidean time. One strategy, proposed in Rothkopf:2009pk and applied for the first time in Rothkopf:2011db to circumvent this problem is to resort to a spectral decomposition of the Wilson loop, which simply amounts to a Fourier transform over a positive definite spectral function

 W□(τ,r)=∫dωe−ωτρ□(ω,r)↔∫dωe−iωtρ□(ω,r)=W□(t,r).

The fact that the time dependence is explicit in the integral kernel tells us that both the real-time and Euclidean time Wilson loop are described by the same spectral function. Hence extracting the values of from imaginary time simulation data will give access to the real-time Wilson loop and in turn to the potential. Inserted into the defining equation (4) the spectrum itself can be related to the values of the potential

 V(r)=limt→∞∫dωωe−iωtρ□(ω,r)/∫dωe−iωtρ□(ω,r). (9)

Unfortunately extracting the spectrum from Euclidean time simulation data is an inherently ill-defined inverse problem, as one seeks to determine the form of a continuous function from a finite and noisy set of individual points. Only by using additional prior information, such as the positive definiteness of the spectrum or smoothness assumptions is it possible to give meaning to this problem, a strategy usually referred to as Bayesian inference. Established implementation of this approach, such as the Maximum Entropy Method Asakawa:2000tr or extended MEM Rothkopf:2012vv have been shown to fail to produce satisfactory results when deployed in the extraction of spectra from Wilson loops Burnier:2013fca , and it needed the development of a new Bayesian reconstruction prescription Burnier:2013nla before quantitatively robust results were obtained. Even with this new method the spectrum can be determined only in a certain range of frequencies, given by the energies resolved on the lattice, which makes a brute force evaluation of (9) impossible. A selection of reconstructed spectra for the two extreme cases of the lowest and highest investigated temperatures around are shown in fig. 1.

This further difficulty can also be overcome from a careful inspection of the spectral structure of the Wilson loop. Just as we did to arrive at eq. (4) let us assume that a potential picture is valid, i.e. the late time evolution of is dominated by a time independent function . As was proven in Burnier:2012az , in this case the Wilson loop spectrum will contain a well defined lowest lying peak of skewed Lorentzian form embedded in a polynomial background

 ρ∝|ImV(r)|cos[Reσ∞(r)]−(ReV(r)−ω)sin[Reσ∞(r)]ImV(r)2+(ReV(r)−ω)2 (10) +c0(r)+c1(r)(ReV(r)−ω)+c2(r)(ReV(r)−ω)2….

The position and width of this peak are related to the real- and imaginary part of the potential respectively. The skewness and background contributions on the other hand are identified with remnants of the physics at scales above the soft-scale.

We can now turn this relation around and use it as a non-perturbative criterion for the applicability of a potential description. If the Wilson loop spectrum does contain a well defined lowest lying peaked feature of skewed Lorentzian type, it will lead to a late-time evolution governed by a well defined potential. As shown in the examples of reconstructed spectra in fig. 1 we do indeed find such well pronounced peaked features and trough a fit establish that they of a skewed Lorentzian form222If the peaks were instead Gaussian , inserting them into the defining formula (9) would lead to a divergent expression, since shows an unphysical linear increase in over time Rothkopf:2011db .. This leads us to conclude that an in-medium potential description is warranted and the values of can be extracted via the application of eq.(10).

In practice, we look for the lowest lying peak in the spectra reconstructed from the lattice QCD observables, fit it around the full width at half maximum with the skewed Lorentzian form (10) and read off the value of and encoded in it. This extraction strategy laid out above has been successfully tested using hard-thermal loop perturbation theory Burnier:2013fca , where both the Wilson loop in Euclidean time, its spectrum and the corresponding potential are known.

In this study we use the values of the heavy quark potential Burnier:2014ssa extracted333Note that the observable we use here to extract the potential is not the Wilson loop but the Wilson line correlator in Coulomb gauge. The reason is that the latter does not contain cusp divergences that make evaluating the Wilson loop very costly on the lattice. We have checked that the result remains gauge invariant by relaxing the Coulomb condition and applying random gauge transformations. from full QCD lattices generated by the HotQCD collaboration Bazavov:2011nk , containing dynamical u,d, and s quarks. The medium quarks on these lattices with spatial extend are described by the asqtad action and are tuned to lie on a line of constant physics with a physical strange quark mass and light u,d quark masses close to their physical values (For more details see Bazavov:2011nk ). In addition to the values of the potential around the phase transition at , which had been extracted in a previous study Burnier:2014ssa , we add here the values from two additional low temperature ensembles close to , which will be used to calibrate the analytic expression for the temperature dependence of the potential in the next subsection. tab. 1 summarizes the relevant simulation parameters for our ensembles and gives the temperatures in relation to the chiral crossover transition temperature on these lattices at  MeV.

Since the finite temperature lattices only contain lattice points in Euclidean time direction, we expect that while a robust determination of spectral peak positions is possible, the accuracy for spectral widths will be insufficient to make reliable statements about . The values obtained for the real-part are given by the colored points in the left panel of fig. 2, the tentative values for from the spectral widths are plotted as lightly colored points in the background of the right panel.

We will postpone the discussion of the temperature dependence of the potential to the end of the next subsection.

### 2.2 Gauss law parametrization of the potential

In order to investigate heavy quarkonium spectra by solving a Schrödinger equation with the proper complex heavy-quark potential, we require an easily evaluable expression for both and at a given temperature. To this end we have recently proposed a field theoretically motivated functional form, that depends on a single temperature dependent parameter, the Debye screening mass Burnier:2015nsa .

The idea behind this approach is as follows. The physics of heavy quarkonium at zero temperature can be described in a potential picture with two distinct characteristics, a Coulombic part at small distances and a linearly rising part at large distances444Note that by allowing the linear term to contribute down to the shortest distances we partially absorb the running of the strong coupling. Hence if we wish to understand the in-medium modification of the inter-quark potential we need to understand how a test charge associated with either one of these two distinct potentials will react to the presence of medium charge carriers in its surroundings. To describe the effect of these charges we will use the hard-thermal loop medium permittivity, in essence making the ansatz of a test particle with a particular electric field configuration being immersed in a bath of weakly interacting quarks and gluons.

The basis for the derivation of the analytic expression for the complex potential lies in the generalized Gauss law derived in Dixit:1989vq

 →∇(→Era+1)=4πqδ(→r). (11)

It is applicable to a point charge with (color) electric field , which covers both the Coulombic , as well as string-like potential relevant for heavy quarkonium. In case of a Coulombic potential we can incorporate the in-medium effects in a straight forward manner by transforming (11) into Fourier space and multiplying the right hand side with a (possibly complex) permittivity, taken here from hard-thermal loop perturbation theory

 ε−1(→p,mD)=p2p2+m2D−iπTpm2D(p2+m2D)2. (12)

Transforming back to coordinate space yields

 −∇2VC(r)+m2DVC(r)=~αs(4πδ(→r)−iTm2Dφ(mDr)), (13)

with

 φ(x)=2∫∞0dpsin(px)pxpp2+1, (14)

which constitutes an integro-differential equation with linear-response character for the in-medium modified Coulomb potential. The strength of the medium effects is clearly governed by the parameter , which we have proposed as a non-perturbative definition of the Debye mass in ref. Burnier:2015nsa . Solving (13) with the appropriate boundary conditions , and reproduces the perturbative potential of Laine et. al Laine:2007qy given in eq.(6).

In the case of a string-like potential, transforming eq.(11) into Fourier space is impractical. Instead we rely on the assumption that the in-medium charge distribution found in the Coulomb case is the same as for the linearly rising potential i.e. a property of the medium and not of the test charge. As was shown in Burnier:2015nsa the strength of the medium effects in the resulting linear-response type equation is governed not by the Debye mass alone but instead by the parameter

 −1r2d2Vs(r)dr2+μ4Vs(r)=σ(4πδ(→r)−iTm2Dφ(mDr)). (15)

Solving for the real part of the medium modified string potential leads to an analytic expression in terms of parabolic cylinder functions

 ReVs(r) =−Γ[14]234√πσμD−12(√2μr)+Γ[14]2Γ[34]σμ. (16)

The imaginary part on the other hand can be written in a closed form using the Wronskian method, which yields

 ImVs(r) = −iσm2DTμ4ψ(μr)=−i~αsTψ(μr), (17)

with

 ψ(x) = D−1/2(√2x)∫x0dyReD−1/2(i√2y)y2φ(ymD/μ) (18) +ReD−1/2(i√2x)∫∞xdyD−1/2(√2y)y2φ(ymD/μ) −D−1/2(0)∫∞0dyD−1/2(√2y)y2φ(ymD/μ).

Note that even though we have used a weak coupling ansatz for the permittivity, which is only appropriate at high temperatures, its insertion into on the Gauss law has produced an expression for the potential with linear-response character. In turn it seems to smoothly connect to zero temperature, as the value of the Debye mass can in principle be set to zero. In particular the real part of the Gauss law derived in-medium potential reduces to the Cornell-potential at . This bodes well for applying the derived expression to fitting the lattice QCD extracted potential values even at temperatures below the formal range of applicability of weak-coupling methods.

### 2.3 Determination of mD from the lattice potential

For a consistent determination of the single temperature dependent parameter , we assume that neither the strong coupling nor the string tension depend on the medium temperature, i.e. all modification emerges from the surrounding light quarks and gluons. Thus we will need to fix the values of , , as well as the arbitrary scale dependent constant shift of at . Since lattice QCD simulations are performed on finite size lattices at a finite lattice spacing, they necessarily operate at a low but finite temperature. In our case we will hence use the newly added data from the two ensembles close to zero temperature at and . We find that the values of the vacuum parameters vary slightly between the two ensembles. Hence a linear interpolation is used at intermediate lattice spacings.555Note that the differences between the values of at different values might be related to an insufficiently precise setting of the scale for the asqtad lattices in Bazavov:2011nk As can be seen from the dark and light blue curves at the top of the right panel of fig. 2 the lattice values for at low temperature indeed show a well pronounced Coulombic and linear behavior that is excellently reproduced by the fit parameters given in tab. 2. I.e. neither the running of the coupling at small distances nor logarithmic corrections at large distances are significant for the regime investigated here.

With the vacuum values set and since the explicit expression for both and derived above only depends on a single parameter, we continue by determining from a fit to the real part of the lattice QCD extracted values alone. The resulting curves are given as solid lines in the right panel of fig. 2, the corresponding Debye masses are collected in tab. 3 and plotted in fig. 3.

We find that tuning allows us to indeed reproduce both the qualitative and quantitative behavior of without problem, if the interpolated constants are used. The strength of the generalized Gauss law Ansatz is that it now allows us to predict the values of the imaginary part of the potential, which was deemed unreliable in a previous study due to the fact that on the finite temperature lattices the Wilson lines were available at only twelve points in Euclidean direction. We are confident in the predictive capabilities of the approach, as it has been able to successfully reproduce in a similar study in quenched QCD Burnier:2015nsa .

When plotting , as implied by via eq. (13) and eq.(17) as solid lines in the right panel of fig. 2 the tentative data (light colored points) shows a surprisingly good agreement with these values at small separation distances fm down to . Only at we find large differences between the extracted and predicted values, which is probably due to the impossibility to resolve the tiny width of the spectral peaks at low temperatures using only data points. Since heavy quarkonium phenomenology in the quark gluon plasma requires knowledge about the potential only until the freezeout boundary is reached, i.e. slightly below the deconfinement transition, the apparent range of applicability of the fit seems to suffice.

### 2.4 Continuum correction for the potential and mD

Several preparations are still in order, since for a meaningful phenomenological investigation, we need to use the continuum values for all parameters entering the potential. In the absence of a true continuum extrapolation of our lattice results we resort to the following strategy.

We fit the vacuum potential parameters, , and the constant term by comparing the energy levels of the corresponding Cornell-potential Hamiltonian to the experimentally measured masses of the bottomonium system, where we expect that finite mass correction are insignificant. More importantly for bottom quarks there exists a well controlled procedure to define their physical mass. A mass often deployed in pNRQCD effective theory calculations is the pole mass Melnikov:2000qh ; Marquard:2007uj ; Gray:1990yh (for bottom GeV) different from the usually quoted mass . Its larger value reflects the fact that the quark mass should be evaluated at the soft scale, relevant for the physics of the bound state and not at the hard scale. A renormalization group flow towards lower energies is hence required, leading to the larger . More precisely, because of the so called renormalon ambiguity in the perturbative evaluation of the Wilson coefficients of the effective theory, one is lead to use the renormalon subtracted scheme Pineda:2001zq . I.e. one subtracts the renormalon ambiguity in the pole mass definition and reshuffles it into the constant part of the potential, which suffers from the same ambiguity. In this way, the ambiguities eventually cancel and we are lead to a quantitatively robust definition of the potential. In general also here one obtains larger values than those resulting from the scheme, such as the established

 mRS′b=4.882±0.041GeV (19)

from ref. Pineda:2001zq , which has been deployed e.g. in Brambilla:2004jw . The errors are estimated to be as large as the highest calculated contribution of order .

This choice of mass allows us to reproduce the PDG values for the four S-wave states , the averaged massed of the three known P-wave triplets as well as the lowest D-wave state to at least three significant digits if the remaining vacuum parameters are set to the following values

 c=−0.1767±0.0210 GeV,~αs=0.5043±0.0298,√σ=0.415±0.015 GeV. (20)

These continuum values are very close to the ones used in conventional quarkonium spectroscopy Quigg:1979vr ; Eichten:1995ch and are also compatible with our lattice data. Note that the effects of string breaking were introduced by hand, flattening off the potential at fm Bali:2005fu . A naive complex scaling analysis confirms that the four lowest S-wave states of the modified Cornell potential Hamiltonian indeed lie below the continuum threshold.

Lattice QCD simulations are carried out in a finite box with a finite lattice spacing, hence discretization effects have to be taken into account. The former is related to the fact that the quarks do not take on their physical masses exactly and thus the chiral crossover temperature on our lattices (MeV) is higher than the continuum value of MeV. In addition the change in lattice spacing leads to slightly different values of the vacuum parameters of the potential in our case, as discussed in sec. 2.3. Hence both effects should be corrected for. Ultimately, we will use the constants of eq.(20) together with the continuum corrected Debye mass

 mphysD(t=T/TphysC)=mD(t)√σ(β)√σcont (21)

where is given in tab. 3 and in (20).

Let us take a closer look at the continuum corrected Debye masses and how they compare to perturbative estimates. In the work of ref. Arnold:1995bh it has been established that the Debye mass can be computed perturbatively only up to the leading order together with the logarithmic correction at next to leading order (NLO). The presence of a magnetic sector in QCD leads to the appearance of truly non-perturbative contributions to at NLO, parametrized in the following by the two terms containing the constants , which need to be determined from numerical simulations

 mD = Tg(μ)√Nc3+Nf6+NcTg(μ)24πlog⎛⎜ ⎜ ⎜ ⎜⎝√Nc3+Nf6g(μ)⎞⎟ ⎟ ⎟ ⎟⎠ (22) +κ1Tg(μ)2+κ2Tg(μ)3.

To compare the temperature dependence of the Debye mass, we fix the non-perturbative constants while keeping constant. For the running of the coupling we utilize the four loop result of ref. Vermaseren:1997fq setting GeV, appropriate for starting the renormalization group flow from a scale where flavors are active.

In the left panel of fig. 3 we show the values of (blue points) together with the fit according to eq. (22). Even though the perturbative part of the formula for is applicable, if at all, at the highest temperatures investigated here, we find that the fit manages to pass through all values of even around the phase boundary. The non-perturbative contributions are non-negligible with and . Nevertheless the perturbation theory motivated fit (by chance) reproduces down to temperatures, slightly below and may hence be used to define a phenomenological, lattice QCD validated, temperature dependence of the Debye screening mass.

## 3 Quarkonium spectra

In sec. 2 we managed to capture the functional form of both the real- and imaginary-part of the static potential by fitting a single temperature dependent parameter , the Debye mass, to extracted on the lattice. We now take the next step and compute from it the in-medium spectral functions of the vector channel bottomonium and charmonium S-wave states at finite mass by solving an appropriate Schrödinger equation. These spectra correspond to resting states with absolute momenta , similar to those usually investigated in direct lattice QCD or lattice NRQCD studies.

The static heavy quark potential is a universal quantity, in the sense that it denotes the lowest order contribution in the non-relativistic expansion for both bottomonium and charmonium physics. Therefore we expect that the vacuum parameters fitted in the bottomonium case do remain the same for the lighter flavor. An important difference between the two cases is that the binding energy of the charmonium ground state is close to . This makes it impossible to set up a perturbative renormalization scheme for the charm mass, similar to the one we used for Bottom. Hence instead of calculating the renormalon subtracted mass, we use the charm mass as fit parameter and tune it to reproduce the masses of the stable S-wave states () and the averaged two P-wave triplets . The resulting best fit value reads

 mPDGfitc=1.472GeV. (23)

Since finite mass effects, in particular radiative corrections can become relevant for charmonium, we expect the agreement between the static potential based masses and the experimental spectrum to be worse than for bottomonium. Indeed for the S-wave states and the 1P triplet only agreement up to the second digit is found when these higher order effects are neglected.

The vacuum parameters determined in sec. 2.4 and (23) constitute the basis from which we embark on the finite temperature study, where all medium effects on the static potential are summarized in the temperature dependence of the Debye mass parameter determined in sec. 2.3. We use the continuum string tension of the bottomonium fit to convert the values of of tab. 3, i.e. the right panel of fig. 3 for use in the continuum Schrödinger equation. Note that we have set up the potential parametrization in sec. 2 such that changing the Debye mass does not affect the overall constant in . I.e. at small enough separation distances, where temperature effects are irrelevant, the values of all agree independently of , which is expected of a correctly renormalized potential and a necessary requirement for a meaningful interpretation of the quarkonium bound state physics at finite . tab. 4 and tab.  5 summarize the vacuum properties of several bottomonium and charmonium states that arise from the potential parameters.

### 3.1 Spectral functions from the Schrödinger equation

All ingredients have been assembled for computing the vector channel spectral functions from the time evolution of the corresponding correlation function governed by the complex in-medium potential666Note again that it is not the Schrödinger equation for the quarkonium wavefunction but for the forward correlator that we are solving here.. In the following we deploy the Fourier space method developed in ref. Burnier:2007qm , which solves

 [~H−i|ImV(r)|]D>(t,r,r′)=i∂tD>(t,r,r′),t>0 (24) [~H+i|ImV(r)|]D>(t,r,r′)=i∂tD>(t,r,r′),t<0 (25)

with

 ~H=2mQ−∇22mQ+Re[V](r)+l(l+1)mQr2 (26)

and the starting condition

 D>(0,r,r′)=−6Ncδ(3)(r−r′). (27)

For the correlator in frequency space we Fourier transform

 ~D(ω,r,r′)≡∫∞−∞dtsiωtD>(t,r,r′), (28)

from which the vector channel spectrum is obtained by taking the limit

 ρV(ω)=limr,r′→012~D(ω,r,r′). (29)

Taking the point splitting to zero requires some effort in the practical implementation as discussed in appendix A of ref. Burnier:2007qm .

In fig. 4 we give overview plots of the computed in-medium S-wave spectral functions for bottomonium and charmonium at several temperatures around the transition temperature. Note that the widths seen here are actual physical widths, related to finite temperature effects, i.e. all bound states reduce to delta peak structures in the limit. By construction, the peaks coincide with the experimental values for the bottomonium states up to 3 digits and fit the two charmonium states below threshold up to few percents 777For , we actually add a small imaginary part in the potential to avoid having exact delta-functions in the spectrum, that would be impossible to visualize graphically..

### 3.2 Properties of the in-medium spectral functions

As expected from the difference in constituent quark mass, the bottomonium states react less severely to the medium at a given temperature compared to charmonium. In general we can observe in this adiabatic setting that the very narrow peaks at low temperature are affected in a hierarchical manner, the highest lying states, which are most loosely bound, begin to broaden and shift first, followed sequentially by the lower lying states. The combination of the effects of screening and scattering encoded in the potential lead to characteristic common changes for both quark flavors. Not only do the bound states broaden but also their mass shifts to lower values. This behavior is found not only for the ground state but for all the higher lying states before they eventually dissolve and become part of the continuum.

To be more quantitative, if a narrow resonance pole lies close to the real frequency axis its spectrum can be described by a simple Breit-Wigner (BW). On the other hand for states that are close to melting, it is necessary to disentangle the remnant bound state signal from the continuum background. For such broad features a skewed Breit-Wigner needs to be considered Taylor , which reads

 ρ(ω≈E)=C(Γ/2)2(Γ/2)2+(ω−E)2+2δ(ω−E)Γ/2(Γ/2)2+(ω−E)2+O(δ2), (30)

where denotes the energy of the resonance, its width and the phase shift.

Using the interpolated form for the Debye mass (22) we perform a temperature scan of the spectrum at every MeV and fit the different peaks with the BW of eq.(30). From these fits we obtain the temperature dependence of the bound state width and mass of different quarkonium states, which are plotted in fig. 5 and fig. 6 respectively. Another quantity particularly relevant to phenomenology is the area under each of the peaks, which we define here as the integrated area of the Breit-Wigner fit function and which is related to the total dilepton emission rate via eq.(8). Its values are given in fig. 7. We find that as temperature increases, the peak area at first remains rather constant, even if the peak becomes wider but eventually and abruptly begins to decrease rapidly towards zero.

We can put these observations in the context of the in-medium modification of the potential shown in fig. 2. While the small distance part of the correctly renormalized is virtually temperature independent, a significant screening of the linear rise at large distances occurs with increasing . It is a peculiarity of the confinement mechanism that with the diminishing remnant of the linear rise; also the threshold to the continuum is lowered. This is reflected in a decrease of the value of (see the gray curve in fig. 5). In turn, the binding energy, defined from the difference between bound state mass and the continuum threshold energy is monotonously lowered. That is to say, the thermal fluctuations of the medium destabilize the bound state. Interestingly, once the threshold moves into the vicinity of a formerly firmly bound quarkonium peak it pushes the spectral feature towards lower frequencies until it eventually disappears.

Let us consider the mass shift observed in fig. 5. At first it might seem counterintuitive, as the expectation for an elementary particle in a medium is exactly the opposite, it will receive a thermal mass. For instance at LO, a single fermion receives a mass correction Donoghue:1984zz ; Seibert:1993aw ; Burnier:2014cna of

 mQ(T)=mQ+δmTQ=mQ+g2T2CF12mQ, (31)

and one might wonder if such a contribution should be added to our potential. As can be seen from (31) it is of higher order in the expansion and hence does not contribute to the static potential. At the next order in the expansion, the potential is corrected by a term of the form Brambilla:2000gk . If the two color charges are widely separated one expects that each of them obtains a thermal mass shift, hence . On the other hand, for one can imagine the neutral bound state will not interact with the plasma so that . The thermal mass shift (31) may hence be considered an upper bound on the thermal mass shift of the bound state. Since it is of higher order it is negligible in our non-relativistic computation, which can be checked using our value of and considering a temperature of MeV. One would obtain MeV for charm and MeV for bottom, insignificant in comparison to the thermal shift observed in fig. 4, which easily reaches MeV (see also fig. 5).

These findings, based on a first principles lattice QCD based complex in-medium potential, are qualitatively similar to what had been observed in potential modeling studies that solve a Schrödinger equation with a potential with imaginary part inserted by hand (see e.g. Petreczky:2010tk ). There only the Coulombic contribution to was used, which leads to more stable behavior than in our case. The presence of the string-like vacuum potential contributes an additional term to , which is of comparable size as the Coulombic contribution to at intermediate temperatures. On the other hand our results differ significantly from the spectra obtained in the T-matrix study of ref. Riek:2010fk . There the authors use purely real model potentials which lead to in-medium states that appear to possess a significantly larger energy than their vacuum counterparts. Updated computations in that framework Liu:2015ypa are expected in the near future. Similar shifts of the resonance peaks to lower energies were also observed in a sum-rule based approach Dominguez:2009mk ; Dominguez:2010mx ; Dominguez:2013fca , but the magnitude of the shift is somewhat weaker there. A more thorough comparison of our results to direct studies in lattice QCD is part of ongoing work.

For MeV one can also compare to the perturbative estimates of Burnier:2008ia and qualitatively good agreement between the spectra is found. The lattice spectra show slightly narrower peaks than the perturbative ones, due to the linearly rising part of the potential, which increase the binding energy. This effect is not captured in perturbation theory. On the other hand, the string effects vanishes quickly at high temperature and in addition the imaginary part of the lattice potential is larger due to string effects. Combined it seems to compensate the stronger binding through an increase in the width of the state. One central benefit of the lattice computation is the possibility to connect the high and low temperature spectra, which is cannot be realized in perturbation theory.

A qualitative difference between the perturbative results and the lattice ones is the way the peak positions change as function of the temperature. In perturbation theory, which contains essentially only the Coulomb term, the bound state peaks move slightly to higher values of frequency Burnier:2007qm as T increases whereas we see here a clear decreases of the peak frequencies.

As the determination of the imaginary part of the potential represents a significant source of uncertainty in this work, we close this section by studying its effect on the spectrum in more detail. To do so, we vary the strength of the imaginary part multiplying it by a -independent factor see fig. 8. The main effect of the imaginary part is to broaden the peak, without changing its position and area, unless the peak is already close to the continuum, as is e.g. the case with in fig. 8. That said, it is obvious that the corresponding states do melt more quickly with a large , as the width reaches the binding energy much more quickly.

## 4 In-medium quarkonium penomenology

### 4.1 Charmonium at freezout

At current heavy-ion colliders, such as RHIC and LHC, with collision energies above GeV, the success of the statistical model of hadronization, to predict the yields of heavy quarkonium, supports the idea that charmonium completely dissolves in the created plasma. In turn essentially all charmonium bound states we observe in experiment would be generated via recombination that takes place at the freeze-out boundary, usually located slightly below the crossover temperature. In such a setting the ratio between the yields for and can be estimated from the difference in area under the corresponding peak structures in the in-medium spectral functions in a straight forward way.

Let us assume that due to the dynamical nature of the expanding fireball the freeze-out happens close to the chiral crossover temperature . When inspecting the spectra at this temperature, as done in fig. 9 we find that there indeed exist two peaked structure related to the in-medium and . I.e. we can compute the amount of in-medium dimuons produced by each of the two states via the product of the spectral function and Bose-Einstein distribution, integrated over the frequency range corresponding to each of the states via eq. (8).

This however is not what is measured in experiment, since the charmonium states do not decay into leptons within the plasma but given their life-time, instead long after the plasma is diluted away. We should thus in principle project the states corresponding to the peaks of the finite spectra onto the states. As no agreed upon method exists to do so, we here assume that the states inside the peaked structures will become real or after freeze out and subsequently decay outside of the plasma. These assumptions are similar to those made in the statistical model Andronic:2009sv and take into account the in-medium modification of the particle states.

In order then to calculate the phenomenologically relevant density of charmonium states at freeze out from the spectral function, we use the following tactic: We estimate first the contribution from the different states to the in-medium dilepton emission rate. From it we can compute the ratio of dileptons produced by and . The ratio of the number densities of and follows when we correct for the probability of the corresponding vacuum states to decay electromagnetically.

In detail, the dilepton emission rate (8) reads

 Rℓ¯ℓ∝∫dp0d3pρ(P)P2nB(p0). (32)

Here we use the fact that to leading order depends only on and leave a more detailed analysis of the momentum dependence for future studies. After performing the change of variable , we obtain

 Rℓ¯ℓ∝∫dωd3pρ(ω)ω2nB(√ω2+p2)ω√ω2+p2. (33)

In this formula, the contribution from the different bound states will come from the corresponding peak area in , see fig. 9. Hence we fit with a skewed Breit-Wigner (30) to distinguish the contribution from the different states . From that fit we obtain the position and width of the spectral feature, i.e. the thermal mass of the particle and its width.

To calculate the integral (33) we checked numerically that it is possible to approximate the Breit-Wigner peak by a delta peak, keeping the area under the curve constant. Performing this, the integral (33) becomes

 RΨnℓ¯ℓ∝A∫dωd3pnB(√M2n+p2)Mn√M2n+p2. (34)

The fit of the charmonium spectrum in fig. 9 then yields

 RΨ′ℓ¯ℓRJ/Ψℓ¯ℓ=0.023±0.004. (35)

To obtain the number density we divide, as discussed, by the electromagnetic decay rate of a vacuum state, which is proportional to the square of the wave function at divided by the square of the mass of the state Bodwin:1994jh

 NΨ′NJ/Ψ∣∣∣T=TC=RΨ′ℓ¯ℓRJ/Ψℓ¯ℓM2Ψ′|ΦJ/Ψ(0)|2M2J/Ψ|ΦΨ′(0)|2=0.052±0.009, (36)

where we used the value at the origin of the zero temperature wave functions. A comparison of our result to those of the statistical model is shown in fig. 10.

While the experimental determination of the to ratio is challenging and currently limited by statistics, both ALICE Adam:2015isa and CMS Khachatryan:2014bva have put forward first measurements of the double ratio . The cancellation of experimental systematic uncertainties makes this quantity the most robust measure of the in-medium relation between and to date. Integrated over centrality CMS found that at the lowest accessible rapidities the ratio lies at , while at larger rapidity it grows to values larger than one.

Here we give a rough estimate of the double ratio by using experimental data obtained for prompt charmonium at TeV by the CMS collaboration. The rapidity averaged cross-section ratios for to denoted by in Chatrchyan:2011kc are furthermore averaged over transverse momentum , as the observed dependence on is rather weak. To extract the number density ratio, the difference in branching fraction into dimuons must also be corrected, after which we obtain

 NΨ′NJ/Ψ∣∣∣√s=7TeV=⟨R⟩pTBR(J/Ψ→μ+μ−)BR(Ψ′→μ+μ−)=0.09±0.015, (37)

This leads to a double ratio of

 NΨ′NJ/Ψ∣∣∣PbPb/NΨ′NJ/Ψ∣∣∣pp=0.58±0.14 (38)

in which the errors in the individual contributions have been naively propagated under the assumption of being independent. Compared to the low rapidity measurements by CMS Khachatryan:2014bva we find good agreement within our still sizable error bars.

### 4.2 Bottomonium at maximum LHC temperatures

Given the small amount of bottom quarks produced in heavy-ion collisions, both at RHIC and LHC, one would naively expect a negligible amount of recombination, a thought which might however need reconsideration Young:2008he ; Young:2009tj . The bottomonium states seen in lead-lead collisions are hence expected to be the ones that survive over the whole plasma evolution. A way to approximate their number is to read off the peak area from the spectrum at the highest temperature reached by the plasma. At LHC, where most experimental results for bottomonium have been obtained, the highest temperature according to ref. Habich:2014jna corresponds to from an analysis of harmonic yields. We hence consider the spectrum at in fig. 11. There, we see that only the ground state survives and is noticeably washed out. Of course not all the bottomonium formed will experience the highest temperature, one e.g. expects states to be also produced at the edge of the plasma, where the temperature is lower.

If we assume the fate of the heavy quarkonium is determined at the maximum temperature reached by the plasma, we can estimate the ratio of bottomonium produced in lead-lead collisions to the amount of Bottom produced in proton-proton collisions, which is given by the ratio of the spectral weights at vs ,

 RPbPb/pp=ΓTmaxCTmaxΓT=0CT=0=0.91±0.02. (39)

Here denotes the width and the amplitude of the corresponding Breit-Wigner type spectral feature. Note that the error bars might be underestimated as the Debye mass had to be extrapolated far above the available temperatures using the HTL fit. The corresponding experimental results from CMS Hu:2011zza are lower with .

### 4.3 Melting temperatures

As we found in fig. 4, the bound states disappear progressively with increasing and we can attempt to define their melting point. In the case of a real potential, the melting temperature is straightforward to define from the disappearance of the purely real binding energy. For a complex potential the situation is more subtle, as the bound state broadens before it disappears. A popular choice is to define the melting of a state the moment its width equals its binding energy Laine:2007qy ; Kakade:2015xua .

Here we scan the spectrum at different temperatures, using the HTL based inter- and extrapolation of the values of the continuum corrected Debye mass . The bound state peak features are determined from a fit based on the skewed Lorentzian of eq.(30). The melting temperatures obtained in this way are given in tab. 6 and exhibit a hierarchical pattern, where only the ground states survive well above 888For charmonium a similar behavior has been anticipated in the in the sequential melting scenario of ref. Karsch:2005nk ..

Our observations are in qualitative agreement with earlier lattice QCD studies of charmonium correlators and spectral functions Ding:2012sp ; Asakawa:2003re ; Datta:2003ww ; Jakovac:2006sf ; Iida:2006mv ; Ohno:2011zc ; Aarts:2007pk ; Borsanyi:2014vka ; Ohno:2014uga . Note that most of these studies were performed in the quenched approximation or at rather large values of the light quark masses and so far no continuum extrapolation was performed. Results of spatial correlation functions and the corresponding spatial screening masses from lattice calculations in the charmoniun sector that include dynamical light quarks further support our findings Bazavov:2014cta ; Karsch:2012na . Relativistic charmonium and bottomonium correlation functions computed in quenched QCD Ohno:2014uga have also shown that in the bottomonium channel thermal modifications set in at larger temperatures, well in the QGP phase, compared to the charmonium channel. The recent determination of bottomonium spectra in full QCD based on non-relativistic QCD Aarts:2014cda ; Kim:2014iga corroborates an ground state survival deep into the QGP. A thorough quantitative comparison of our lattice potential based spectra to direct lattice QCD computations will be part of a future study and is work in progress.

A more precise determination of the melting temperatures will require lattice simulations at higher temperatures, as for now we can only use an extrapolation of the HTL fit for , which becomes increasingly unreliable at high temperatures beyond .

## 5 Conclusion

Heavy quarkonium is a vital probe to uncover the physics of the quark-gluon plasma created in heavy-ion collisions. Their non-relativistic nature opens up the possibility to describe their in-medium behavior by an effective potential entering a Schrödinger equation. We reported here the first quantitative phenomenological investigation of bottomonium and charmonium S-wave spectra at finite , based on the proper complex in-medium static potential extracted from lattice QCD.

To make possible the use of the lattice potential in a phenomenological application we deployed the generalized Gauss law ansatz, which provides analytic expressions for and that depend on a single temperature dependent parameter , the Debye mass. Tuning it was possible to reproduce the lattice values of and to validate the up to now tentative values of at all separation distances and temperatures investigated. After correcting for lattice discretization artifacts in the Debye mass we fitted its temperature dependence successfully with a HTL based expression.

To compute spectral functions for phenomenological inspection, we determined the vacuum parameters of the static potential in the continuum from a fit to the experimentally known S-, P- and D-wave bottomonium states. Finite effects are incorporated though the continuum corrected Debye mass, extracted from the real-part of the in-medium lattice potential, leading to correctly renormalized finite temperature expressions for and

The central findings in our adiabatic setup are a clear pattern of sequential melting of both bottomonium and charmonium with respect to their vacuum binding energies. The interplay between Debye screening and scattering with medium partons leads to characteristic shifts of the in-medium bound states to lower masses before they dissolve into the continuum. This effect is opposite to the relatively small gain in thermal mass.

The availability of in-medium spectra with physical widths allowed us to estimate the to ratio at the crossover transition temperature, relevant in the context of the proposed statistical hadronization scenario. Our approach to compute the ratio of in-medium dimuon emission corrected by the emission rate of the corresponding vacuum states yields a value slightly larger than that of the statistical model but still much smaller than the experimental data in collisions. A rough estimate of the suppression of the state in this adiabatic scenario on the other hand gave a value larger than experimentally observed.

The static potential is only the starting point for a quantitatively reliable investigation of heavy quarkonium. In particular for the lighter flavor charmonium the question of finite mass effects in the potential should be addressed and efforts should be directed at determining the first correction at finite temperature from the lattice. On the other hand direct lattice QCD studies, be it in the relativistic formulation or using the effective field theory NRQCD including velocity correction, will be essential to crosscheck the validity of the potential description. Comparisons with the data here can be made both on the level of spectra, as well as the corresponding Euclidean correlators. The necessity to solve an inherently ill-defined problem in order to extract the spectral functions from lattice QCD current-current correlators directly remains a significant challenge. In order to surmount it both methodological progress in Bayesian and non-Bayesian spectral reconstruction approaches, as well as computational efforts to generate lattices with more Euclidean time steps will be required.

The lattice in-medium potential can also be used to setup a description of the real-time evolution of the heavy quarkonium wave function at finite temperature in the context of open quantum systems. An approach based on a stochastic potential Akamatsu:2011se appears promising, where a purely real potential, given by is randomly perturbed at each time step by noise of strength related to . If developed further it promises to provide vital and phenomenologically relevant information on e.g. the density matrix of states for the quarkonium system, which goes beyond what is accessible from within the spectral functions computed here.

Directions for future work include the determination of the complex in-medium potential on lattices at higher temperature in order to avoid the extrapolations beyond required e.g. in the determination the melting temperature in sec. 4.2. Another aspect of interest is the momentum dependence of the spectra which is also an active area of research in direct lattice QCD studies.

## Acknowledgments

The authors thank G. Aarts, A. Andronic, N. Brambilla, M. Laine, P. Petreczky and B. Keren-Zur for stimulating discussions. We furthermore thank W. Soeldner and the Bielefeld-BNL-CCNU Collaboration for providing the asqtad gaugefield configurations used in this study. Part of the calculations were performed on the Bielefeld GPU cluster and on the ITP in-house cluster in Heidelberg. YB is supported by SNF grant PZ00P2-142524.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters