Quarkonium at finite temperature: Towards realistic phenomenology from first principles

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 ${\rm Re}[V]$ with a single temperature dependent parameter $m_D$, the Debye screening mass, and confirm the up to now tentative values of ${\rm Im}[V]$. The obtained analytic expression for the complex potential allows us to compute quarkonium spectral functions by solving an appropriate Schr\"odinger 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 (${\rm Re}[V]$) and scattering (${\rm Im}[V]$) 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.


Introduction
Heavy quarkonium, the bound states of a heavy quark and anti-quark (cc or bb), 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 [1]. 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 [2,3]).
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 T = 0 in [4] and investigated in the non-perturbative regime in [5,6], with both cases being reviewed in detail in [7]. The generalization to finite temperature in a perturbative setting followed in [8].
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 v, to setup a simplified description of the bound state evolution at soft energy scales (E soft ∼ m Q v) 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 (E hard ∼ 2m Q ) 1 . And indeed the heavy quark rest mass (in this work we use m c 1.472GeV and m b = 4.882GeV) and the characteristic scale of quantum fluctuations in QCD, usually denoted by Λ QCD 200MeV 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 QQ can be described by the thermally averaged unequal-time point-split meson-meson correlator with the meson operator defined as M (x, y, t) =Q(x, t)γ µ U [x, y]Q(y, t). γ µ refers to the Dirac matrices and U [x, y] denotes a straight Wilson line connecting (x, t) and (y, t). The relative coordinate enters as r = x − y.
In the limit of infinite quark mass, where the static constituents are truly test particles, (1.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 r = |r| becomes an external parameter of the theory and their evolution traces out a rectangle over time. Formally it has been shown that eq.(1.1) reduces to the rectangular real-time Wilson loop [7] W (r, t) = Tr exp − ig dx µ A a µ T a . (1.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 [10] i∂ t W (r, t) = Φ(r, t)W (r, t) (1. 3) with Φ(r, t) being a time-and space dependent complex function.
In case that a potential picture for the static QQ system is applicable, the function Φ(r, t) has to asymptote to a time independent function V (r). In turn it will dominate the evolution of W (r, t) at times much larger than the intrinsic scales of e.g. the gluons and light quarks in the medium. Formally we may then write . (1.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 T 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. [1,5,6]), 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. [7,11]). In a nonperturbative setting, if non-potential effects become sizable, the late time evolution of the Wilson loop cannot be described by a simple time-independent potential V (r). 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 T = 0. 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 [12], readily calculable in lattice QCD. In particular two quantities have gained popularity as model potentials, the color singlet free energies in Coulomb gauge defined from the correlator of Polyakov loops Ω(r) and a derived quantity the color singlet internal energies U (1) = F (1) − T S. 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 [15,16]. 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. [8,10] in 2007. In their study the real-time definition (1.4) was evaluated in a resummed perturbative framework, called the hardthermal 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 Here we absorb a factor C F in the definition of the coupling constantα s = g 2 C F 4π 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 colorsingle free energies in hard-thermal loop resummed perturbation theory in [15,16], 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 [10,17] with the heavy quarks (Landau damping) related to Im[V ] 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 [8] and in turn to the dissociation of the QQ system into gluons [18].
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, Q andQ 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. [19]) 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.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 T = 0 broaden and shift as screening and scattering modifies the QQ state with increasing temperature. Choosing a popular criterion of melting temperature [10,60] 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 [13] dR where Q q denotes the electric charge of the heavy quark considered in units of e, the four momentum is P = (p 0 , p) and the finite mass of the leptons has been neglected (for a more detailed discussion of the above formula see ref. [14]). 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. [30], which provides an analytic formula for Re [V] and Im[V] depending on a single temperature dependent parameter m D the Debye mass. Tuning m D we are able to reproduce the lattice values of Re[V] and confirm the yet tentative values for Im [V]. These values for m D 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 Swave channel and investigate their modification with increasing temperature. Besides giving the melting temperatures for different states, we will determine the Ψ to J/Ψ ratio at the chiral crossover and give a rough estimate for the suppression of bottomonium in a heavy-ion collision compared to p + p. We close with a conclusion in sec. 5.

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 [20]. This calls for an evaluation of the potential definition (1.4) in lattice QCD, which at first seems unfeasible, since direct access to real-time quantities, such as the Wilson loop (1.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.

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 [21] and applied for the first time in [22] 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 ρ 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 (1.4) the spectrum itself can be related to the values of the potential 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 [23] or extended MEM [24] have been shown to fail to produce satisfactory results when deployed in the extraction of spectra from Wilson loops [25], and it needed the development of a new Bayesian reconstruction prescription [26] 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 (2.1) impossible. A selection of reconstructed spectra for the two extreme cases of the lowest and highest investigated temperatures around T C 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. (1.4) let us assume that a potential picture is valid, i.e. the late time evolution of W (t, r) is dominated by  2) up to quadratic polynomial order, we fit the points in the dark red region close to the maximum in each case. The resulting fit functions (dashed orange) reproduce the shape of the spectra (blue circles) quantitatively far beyond the actual fitting region. Note that this indicates that the peak is a Lorentzian and not e.g. a Gaussian a time independent function lim t→∞ Φ(t, r) = V (r). As was proven in [27], in this case the Wilson loop spectrum will contain a well defined lowest lying peak of skewed Lorentzian form embedded in a polynomial background 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 form 2 . This leads us to conclude that an in-medium potential description is warranted and the values of V (r) can be extracted via the application of eq.(2.2).
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 (2.2) and read off the value of Re[V ] and Im[V ] encoded in it. This extraction strategy laid out above has been successfully tested using hard-thermal loop perturbation theory [25], 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 [28] extracted 3 from full QCD lattices generated by the HotQCD collaboration [29], containing dynamical u,d, and s quarks. The medium quarks on these lattices with spatial extend N s = 48 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 m s /m u,d = 20 (For more details see [29]). In addition to the values of the potential around the phase transition at T /T C = 0.68 − 1.66, which had been extracted in a previous study [28], we add here the values from two additional low temperature ensembles close to T ∼ 0, 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 T C = 172.5 MeV.
Since the finite temperature lattices only contain N τ = 12 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 Im [V ]. The values obtained for the real-part are given by the colored points in the left panel of fig. 2, the tentative values for Im[V ] from the spectral shows an unphysical linear increase in Im[V] over time [22]. 3 Note 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. 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.

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 Re[V ] and Im[V ] 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 [30].
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 distances 4 . 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 [31] ∇ E r a+1 = 4π q δ( r). (2. 3) It is applicable to a point charge with (color) electric field E = qr a−1r , which covers both the Coulombic a = −1, q =α s , [α s ] = 1, as well as string-like potential a = 1, q = σ, [σ] = GeV 2 relevant for heavy quarkonium. In case of a Coulombic potential we can incorporate the in-medium effects in a straight forward manner by transforming (2.3) into Fourier space and multiplying the right hand side with a (possibly complex) permittivity, taken here from hard-thermal loop perturbation theory (2.4) Transforming back to coordinate space yields 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 m D , which we have proposed as a non-perturbative definition of the Debye mass in ref. [30]. Solving (2.5) with the appropriate boundary conditions ReV C (r)| r=∞ = 0, ImV C (r)| r=0 = 0 and ∂ r ImV C (r)| r=∞ = 0 reproduces the perturbative potential of Laine et. al [10] given in eq.(1.6).
In the case of a string-like potential, transforming eq.(2.3) 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 [30] 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 Solving for the real part of the medium modified string potential leads to an analytic expression in terms of parabolic cylinder functions D ν (x) The imaginary part on the other hand can be written in a closed form using the Wronskian method, which yields  with 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 T = 0 Cornell-potential at m D = 0. 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.

Determination of m D from the lattice potential
For a consistent determination of the single temperature dependent parameter m D , 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α s , σ, as well as the arbitrary scale dependent constant shift of Re[V ] at T = 0. 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 β 1 = 6.9 and β 2 = 7.48. 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. 5 As can be seen from the dark and light blue curves at the top of the right panel of fig. 2  With the vacuum values set and since the explicit expression for both Re[V ] and Im[V ] derived above only depends on a single parameter, we continue by determining m D 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 m D allows us to indeed reproduce both the qualitative and quantitative behavior of Re[V ] without problem, if the interpolated T = 0 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 β 6.8 6 Table 3. Debye masses extracted from the isotropic HotQCD 48 3 × 12 lattices with asqtad action. For use in phenomenology, a continuum corrected m D may be obtained from the ratio m D / σ(β) shown here, through a multiplication with the continuum value of σ. the predictive capabilities of the approach, as it has been able to successfully reproduce Im[V ] in a similar study in quenched QCD [30].
When plotting Im[V ], as implied by m D (T ) via eq. (2.5) and eq.(2.9) 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 r < 0.75fm down to T = 0.95T C . Only at T = 0.86T C 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 N τ = 12 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.

Continuum correction for the potential and m D
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,α s , σ and the constant term c 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 [32][33][34]  . 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 [35]. 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 M S scheme, such as the established m RS b = 4.882 ± 0.041 GeV (2.11) from ref. [35], which has been deployed e.g. in [7]. These continuum values are very close to the ones used in conventional quarkonium spectroscopy [2,3] and are also compatible with our lattice data. Note that the effects of string breaking were introduced by hand, flattening off the potential at r sb = 1.25fm [36]. 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 (T lat C = 172.5MeV) is higher than the continuum value of T cont C = 155 ± 9MeV. 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.(2.12) together with the continuum corrected Debye mass m phys is given in tab. 3 and σ cont in (2.12).
Let us take a closer look at the continuum corrected Debye masses and how they compare to perturbative estimates. In the work of ref. [37] 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 m D at NLO, parametrized in the following by the two terms containing the constants κ 1 , κ 2 , which need to be determined from numerical simulations (2.14) To compare the temperature dependence of the Debye mass, we fix the non-perturbative constants κ 1 , κ 2 while keeping µ = 2πT constant. For the running of the coupling g(µ) we utilize the four loop result of ref. [38] setting Λ QCD = 0.2145GeV, appropriate for starting the renormalization group flow from a scale where N f = 5 flavors are active.
In the left panel of fig. 3 we show the values of m D /T (blue points) together with the fit according to eq. (2.14). Even though the perturbative part of the formula for m D is applicable, if at all, at the highest temperatures investigated here, we find that the fit manages to pass through all values of m D even around the phase boundary. The non-perturbative contributions are non-negligible with κ 1 = 0.84 ± 0.10 and κ 2 = −0.40 ± 0.03. Nevertheless the perturbation theory motivated fit (by chance) reproduces m D down to temperatures, slightly below T c and may hence be used to define a phenomenological, lattice QCD validated, temperature dependence of the Debye screening mass.

Quarkonium spectra
In sec. 2 we managed to capture the functional form of both the real-and imaginarypart of the static potential by fitting a single temperature dependent parameter m D , the Debye mass, to Re[V ] 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 p = 0, 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 Λ QCD . 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 (J/Ψ, Ψ ) and the averaged two P-wave triplets χ c (1P ), χ c (2P ). The resulting best fit value reads m PDG fit c = 1.472 GeV. (3.1) 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 T = 0 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.    Table 5. The masses, mean radii and distances to the DD threshold for charmonium S-wave and P-wave states at T=0.

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 D > (t, r, r ) governed by the complex in-medium potential 6 . In the following we deploy the Fourier space method developed in ref. [39], which solves and the starting condition For the correlator in frequency space we Fourier transform (3.6) from which the vector channel spectrum is obtained by taking the limit Taking the point splitting to zero requires some effort in the practical implementation as discussed in appendix A of ref. [39].
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 T = 0 limit. By construction, the T = 0 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 7 .

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 Re[V] and scattering Im[V] 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 [50], which reads The charmonium in-medium spectral function based on the same complex potential as for bottomonium. In the absence of a well defined renormalon subtracted charm mass, its value is fitted to reproduce the known charmonium bound states at T = 0 giving m c = 1.472GeV. Also charmonium shows a sequential melting pattern, as we find that Ψ(2S) disappear around T C while J/Ψ close to 1.41T C .
where E denotes the energy of the resonance, Γ its width and δ the phase shift.
Using the interpolated form for the Debye mass (2.14) we perform a temperature scan of the spectrum at every δ T = 3 MeV and fit the different peaks with the BW of eq.(3.8). 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 A = πCΓ 2 and which is related to the total dilepton emission rate via eq.(1.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 Re[V] is virtually temperature independent, a significant screening of the linear rise at large distances occurs with increasing T . It is a peculiarity of the confinement mechanism that with the diminishing remnant of the linear rise; also the threshold to the continuum E cont is lowered. This is reflected in a decrease of the value of E cont = Re[V(r → ∞)] (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 [40][41][42] of and one might wonder if such a contribution should be added to our potential. As can be seen from (3.9) it is of higher order in the 1/m Q expansion and hence does not contribute to the static potential. At the next order in the 1/m Q expansion, the potential V (r) is corrected by a term of the form V 1 (r)/m Q [5]. If the two color charges are widely separated one expects that each of them obtains a thermal mass shift, hence lim r→∞ V 1 (r)/m Q = 2δm T Q . On the other hand, for r → 0 one can imagine the neutral bound state will not interact with the plasma so that V 1 (r) = 0. The thermal mass shift (3.9) 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α S = g 2 C F /4π and considering a temperature of T ∼ 200MeV. One would obtain δm T Q = 7MeV for charm and δm T Q = 2MeV for bottom, insignificant in comparison to the thermal shift observed in fig. 4, which easily reaches 100 − 300MeV (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. [43]). There only the Coulombic contribution to Im[V ] was used, which leads to more stable behavior than in our case. The presence of the stringlike vacuum potential contributes an additional term to Im[V ], which is of comparable size as the Coulombic contribution to Im[V ] at intermediate temperatures. On the other hand our results differ significantly from the spectra obtained in the T-matrix study of ref. [44]. 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 [45] are expected in the near future. Similar shifts of the resonance peaks to lower energies were also observed in a sum-rule based approach [46][47][48], 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 T ∼ 250MeV one can also compare to the perturbative estimates of [49] 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 [39] as T increases whereas we see here a clear decreases of the peak frequencies. Ρ Ω Ω 2 Figure 8. Charmonium (left) and Bottomonium (right) spectrum at T = 1.19T C with manually changed values for the imaginary part of the potential (0.1 ImV (black), ImV (red) and 2 ImV (blue)). While the peak position hardly changes unless the peak is close to melting, the width significantly depends on the strength of ImV).
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 r-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 Υ(3S) in fig. 8. That said, it is obvious that the corresponding states do melt more quickly with a large Im[V ], as the width reaches the binding energy much more quickly.

.1 Charmonium at freezout
At current heavy-ion colliders, such as RHIC and LHC, with collision energies above √ s AA > 200GeV, 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 J/Ψ and Ψ(2S) 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 freezeout happens close to the chiral crossover temperature T C . 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 J/Ψ 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. (1.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 T spectra onto the T = 0 states. As no agreed upon method exists to do so, we here assume that the states inside the peaked structures will become real J/Ψ or Ψ after freeze out and subsequently decay outside of the plasma. These assumptions are similar to those made in the statistical model [51] 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 J/Ψ. The ratio of the number densities of Ψ and J/Ψ follows when we correct for the probability of the corresponding vacuum states to decay electromagnetically.
In detail, the dilepton emission rate (1.8) reads R ¯ ∝ dp 0 d 3 p ρ(P ) P 2 n B (p 0 ). (4.1) Here we use the fact that to leading order ρ depends only on P 2 = p 2 0 − p 2 and leave a more detailed analysis of the momentum dependence for future studies. After performing the change of variable ω = p 2 0 − p 2 , we obtain In this formula, the contribution from the different bound states will come from the corresponding peak area in ρ(ω)/ω 2 , see fig. 9. Hence we fit ρ(ω)/ω 2 with a skewed Breit-Wigner (3.8) to distinguish the contribution from the different states n. From that fit we obtain the position and width of the spectral feature, i.e. the thermal mass M n of the particle and its width.
To calculate the integral (4.2) we checked numerically that it is possible to approximate the Breit-Wigner peak by a delta peak, keeping the area A under the curve constant. Performing this, the integral (4.2) becomes The fit of the charmonium spectrum in fig. 9 then yields 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 r = 0 divided by the square of the mass of the state [52] 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 J/Ψ ratio is challenging and currently limited by statistics, both ALICE [53] and CMS [54] have put forward first measurements of the double ratio  Figure 10. Comparison of predicted yields for Ψ and J/Ψ from this work (orange) compared to predictions from the statistical model (green) and several experimental measurements. The plot structure and experimental data has been adapted from ref. [51].
systematic uncertainties makes this quantity the most robust measure of the in-medium relation between Ψ and J/Ψ to date. Integrated over centrality CMS found that at the lowest accessible rapidities the ratio lies at 0.45 ± 0.13(stat) ± 0.07(syst) , while at larger rapidity 1.6 < |y| < 2.4 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 √ s = 7TeV by the CMS collaboration. The rapidity averaged cross-section ratios for Ψ to J/Ψ denoted by R in [55] are furthermore averaged over transverse momentum R p T = 0.0378 ± 0.006, as the observed dependence on p T is rather weak. To extract the number density ratio, the difference in branching fraction into dimuons must also be corrected, after which we obtain This leads to a double ratio of the assumption of being independent. Compared to the low rapidity measurements by CMS [54] we find good agreement within our still sizable error bars.

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 [56,57]. 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. [58] corresponds to T 3T C from an analysis of harmonic yields. We hence consider the spectrum at T max = 3T C 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 T max = 3T C vs T = 0, Here Γ denotes the width and C 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 [59] are lower with R PbPb/pp = 0.6 ± 0.2.

Melting temperatures
As we found in fig. 4, the bound states disappear progressively with increasing T 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 [10,60].
Here we scan the spectrum at different temperatures, using the HTL based interand extrapolation of the values of the continuum corrected Debye mass m D . The bound state peak features are determined from a fit based on the skewed Lorentzian of eq.(3.8). 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 T C 8 . Our observations are in qualitative agreement with earlier lattice QCD studies of charmonium correlators and spectral functions [62][63][64][65][66][67][68][69][70]. 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 [71,72]. Relativistic charmonium and bottomonium correlation functions computed in quenched QCD [70] 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 [73,74] corroborates an Υ(1S) 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 Υ(1S) melting temperatures will require lattice simulations at higher temperatures, as for now we can only use an extrapolation of the HTL fit for m D , which becomes increasingly unreliable at high temperatures beyond T > 1.66T C .   Table 6. Melting temperatures T melt of the different bound states defined by the point, at which the width of the state equates its binding energy calculated from Re[V] alone. The error on the determination of T melt takes into account the possible variation of m D as shown in fig. 3. Note that because of the lack of sufficient number of lattice ensembles below T = 0.95T C and the breakdown of the generalized Gauss law ansatz, we only give upper limits in this regime.

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 T , based on the proper complex in-medium static potential extracted from N f = 2 + 1 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 Re [V] and Im[V] that depend on a single temperature dependent parameter m D , the Debye mass. Tuning m D it was possible to reproduce the lattice values of Re [V] and to validate the up to now tentative values of Im[V] 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 T 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 Re [V] and Im [V] 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 J/Ψ 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 p + p collisions. A rough estimate of the suppression of the Υ(1S) state in this adiabatic scenario on the other hand gave a value 50% 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 V 1 (r)/m Q 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 [19] appears promising, where a purely real potential, given by Re[V] is randomly perturbed at each time step by noise of strength related to Im[V]. 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 T = 1.66T C required e.g. in the determination the Υ(1S) 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.