Improved Gauss law model and in-medium heavy quarkonium at finite density and velocity

We explore the in-medium properties of heavy-quarkonium states at finite baryo-chemical potential and finite transverse momentum based on a modern complex valued potential model. Our starting point is a novel, rigorous derivation of the generalized Gauss law for in-medium quarkonium, combining the non-perturbative physics of the vacuum bound state with a weak coupling description of the medium degrees of freedom. Its relation to previous models in the literature is discussed. We show that our approach is able to reproduce the complex lattice QCD heavy quark potential even in the non-perturbative regime, using a single temperature dependent parameter, the Debye mass mD. After vetting the Gauss law potential with state-of-the-art lattice QCD data, we extend it to the regime of finite baryon density and finite velocity, currently inaccessible to first principles simulations. In-medium spectral functions computed from the Gauss law potential are subsequently used to estimate the ψ′/J/ψ ratio in heavy-ion collisions at different beam energies and transverse momenta. We find qualitative agreement with the predictions from the statistical model of hadronization for the √ sNN dependence and a mild dependence on the transverse momentum. Corresponding author. ar X iv :1 90 6. 00 03 5v 2 [ he pph ] 1 3 Ju n 20 19


Introduction
Heavy quarkonium, the bound states of a charm or bottom quark and its anti-quark (cc or bb), have matured into a high precision tool for the study of strongly interacting matter under extreme conditions in the context of relativistic heavy-ion collisions [1,2]. On the experimental side, the predominant decay of heavy quarkonia into dileptons makes them well-controlled observables, probing different stages of the quark gluon plasma (QGP) created in such collisions. The now iconic measurements of the dimuon spectrum of bottomonium by the CMS collaboration [3], lend themselves to a phenomenological interpretation of the bb pair traversing the medium as a test particle (see for example [4]), sampling the whole history of the QCD medium. On the other hand the more recent measurement of a finite elliptic flow of the J/ψ particle by the ALICE collaboration [5] indicates at least a partial equilibration of the charm quarks with their surrounding. The inevitable loss of memory about the initial conditions accompanying equilibration hence positions the lighter flavor as a probe of the late stages of the collision.
From a theory standpoint the heavy mass of the constituent quarks opens the door to the powerful effective field theory (EFT) framework [6] that allow us to simplify the description of their (non-)equilibrium behavior. These techniques have led to progress both in direct lattice QCD studies of equilibrated quarkonium as well as in formulating real-time descriptions of their non-equilibrium evolution. The foundation of the EFT strategy is the presence of the natural separation of scales m Q m Q v m Q v 2 , with m Q the heavy quark mass and v its typical velocity. These are denoted the hard (rest-energy), soft (momentum exchange) and ultra-soft (binding energy) scales respectively. In addition to these three scales, there exists the characteristic scale of quantum fluctuations Λ QCD and of thermal fluctuations T .
The construction of an EFT involves first choosing a cut-off energy above which the physics is not treated explicitly, before identifying the relevant degrees of freedom in the lower energy scale and treating these explicitly by writing down the most general Lagrangian compatible with the corresponding symmetries. Each term receives a complex-valued Wilson coefficient which needs to be determined by a matching procedure; that is, we compute a correlation function in the full microscopic theory and a correlator with same physics content in the EFT and require that they agree below a certain scale. By this integrating out of higher scales, the results of the microscopic theory can be reproduced in the EFT as long as we stay within its range of validity at low energies.
For heavy quarkonium this program has been implemented by integrating out the hard scale ∼ m Q from the full QCD Lagrangian to give Non-Relativistic QCD (NRQCD), a theory of non-relativistic Pauli spinor fields. This can be achieved in a fully non-perturbative manner. In a further step, integrating out the soft scale ∼ m Q v results in Potential Non-Relativistic QCD (pNRQCD) [7], where the potential governing the quarkonium dynamics becomes one of the Wilson coefficients determined via matching. While the perturbative derivation of pNRQCD has been successfully completed, its non-perturbative definition is still an active field of research. In this work we will utilize the fact that the in-medium potential between two static quarks can indeed be systematically derived from QCD in a non-perturbative fashion via the language of EFTs. 1 Much progress has been made in understanding the properties of equilibrated heavy quarkonium in a static medium from first principles by using lattice QCD computations. To directly study the in-medium modification of charmonium it is nowadays possible to deploy fully relativistic formulations of the heavy quarks (for some recent works see [10][11][12]). However, realistic simulations of bottomonium currently deploy lattice regularized versions of NRQCD [13,14]. In-medium correlation functions computed from first principles have revealed the presence of statistically significant in-medium modifications consistent with our intuition: the hotter the medium, the stronger the modification, and the more weakly bound the quarkonium state is in vacuum, the easier it is influenced by the medium.
Going beyond statements of overall in-medium modification remains difficult in lattice simulations, as it requires one to extract the spectral functions of in-medium quarkonium from Euclidean time correlators. This amounts to an exponentially hard ill-posed inverse problem, which despite concerted efforts of the community so far has only revealed insight into the properties of the ground state in-medium properties of different quarkonium channels. The most recent study of bottomonium and charmonium from lattice NRQCD has elucidated the change in mass of the ground state in-medium, concluding that quarkonium becomes lighter as temperature increases [14].
What all of these studies are still missing is access to the remnants of excited states in the medium, as well as to the threshold, which plays an important role in understanding the stability of the in-mediums states. Indeed the in-medium binding energy is defined from the distance of the in-medium quarkonium spectral peak to the onset of the threshold. Currently this information is only accessible in potential based computations, where a Schrödinger equation for the spectral functions is solved in the presence of a non-perturbatively defined in-medium potential.
Important progress in this regard has been made using an EFT based definition of the in-medium potential between two static quarks based on the real-time evolution of the QCD Wilson loop, . (1.1) Evaluating this expression in hard-thermal loop (HTL) resummed perturbation theory revealed [6,15,16] that in general the proper potential is a complex quantity with a real part exhibiting Debye screening and an imaginary part growing monotonously with temperature. The physical processes contributing to ImV differ according to the separation of scales present [17]; the scattering of medium partons with the gluons mediating the heavy quark interaction, so called Landau damping, as well as the gluon induced transition from a color singlet to octet may both contribute. At high temperature, Landau damping dominates and the potential reads Note that this potential describes the real-time evolution of an unequal time quarkonium correlation function and not of the wavefunction itself. Therefore the presence of an imaginary part is not directly related to the disappearance of the heavy quarks from the system but instead encodes the decoherence of the system from its initial state as it evolves over time in the thermal medium [18,19]. It is important to keep in mind that when we solve a Schrödinger equation with this potential, the resulting correlation function can be used to straightforwardly compute the in-medium quarkonium spectral function. The question of how to relate this complex-valued potential to the real-time evolution of the wavefunction is an active field of research, and recent progress has been made by considering the concept of open-quantum-systems [20][21][22][23][24].
The real-time definition Eq. (1.1) is not directly amenable to an evaluation in lattice simulations since they are carried out in an artificial Euclidean (imaginary) time. Instead a strategy has been developed to use Bayesian inference to extract the values of the potential non-perturbatively on the lattice [25][26][27]. Many studies have confirmed (see for example [28]) that at low temperatures the potential is well described by the Cornell form whereα s = C F g 2 /4π is the strong coupling constant re-scaled to match phenomenology literature, σ is the string-tension and c an additive constant. Eq. (1.4) already captures the two most important features of QCD: asymptotic freedom via the running coupling at small distances and confinement via the non-perturbative linear rise. At finite temperature lattice QCD tells us that not only the real part weakens gradually as one moves into the deconfined phase but that above the pseudo-critical temperature indeed a finite imaginary part is present [28,29]. To put these numerical results to use in computations of quarkonium spectral functions an efficient analytic parametrization of the complex valued potential is needed. Deploying it in a Schrödinger equation gives access to spectral functions from which we may learn about physically relevant properties, such as in-medium masses or decays widths.
To this end, in [30] one of the authors proposed a simple model of a non-perturbative vacuum bound state immersed in a weakly coupled medium (see also [31,32]). The former is described by the Cornell potential, the latter by an in-medium permittivity evaluated in HTL perturbation theory. In that model the effects of the medium on the vacuum potential are incorporated by application of the generalized Gauss law [33]. Once the vacuum parameters of the Cornell potential are chosen, the model provides a prediction for the full in-medium values of ReV and ImV based on a single temperature dependent parameter, identified as the Debye mass m D . While it is not obvious that such an ansatz can accommodate the physics of the heavy quark potential, especially in the non-perturbative regime close to the crossover temperature, it has been shown that tuning of m D reproduces the lattice QCD values of ReV and the tentative values of ImV quite well. In turn the Gauss law model has been used to study the in-medium properties of quarkonium in a thermal medium, as well as to provide an estimate for the ψ /J/ψ ratio at very high energy heavy-ion collisions at mid-rapidity and zero transverse momentum.
However, there exists two shortcoming of the Gauss law model of [30], one technical and one phenomenological, which limit its utility in exploring heavy quarkonium in heavy-ion collisions. The technical one is related to the fact that in order to derive the in-medium modification of the string part of the Cornell potential, the previous study introduced an ad-hoc assumption about the functional form of the in-medium permittivity in coordinate space. On the phenomenological side, the model did not incorporate the effects of finite baryo-chemical potential or transverse momentum, both of relevance to compare to actual data from heavy-ion collisions.
The present study sets out to overcome both of these issues. As a first step we put forward a novel and improved Gauss law model, based on a more rigorous derivation, which does not rely on any ad-hoc ingredients. Taking into account string breaking, i.e. the fact that the vacuum Cornell potential does not rise indefinitely, we are able to derive well defined expressions for ReV and ImV . Their functional form turns out to be simpler than the one obtained in [30]. Again, ReV and ImV only depend on a single temperature dependent parameter m D . They furthermore consistently reduce to the HTL result at large values of the Debye mass parameter (high temperature) and to the Cornell potential at vanishing m D (vacuum). Using state-of-the-art lattice QCD results for ReV and ImV we show that the new model reproduces the non-perturbative values excellently by an appropriate selection of m D . The relation of this new Gauss law model to previous model potentials in the literature is also discussed.
The second improvement is related to extending the model to settings not accessible to first principle lattice QCD simulations, i.e. quarkonium in a medium at finite baryochemical potential, as well as quarkonium traversing the QGP with a finite velocity. The latter is implemented via the hot-wind scenario, where the medium moves with a finite relative velocity with respect to the quarkonium [34]. The extended model will be used to compute the in-medium spectral functions of both charmonium and bottomonium states and investigate their in-medium properties as well as their melting. Modeling the finite baryo-chemical potential regime of the QCD phase diagram allows us to estimate the ψ /J/ψ ratio in heavy-ion collisions at lower beam energies, relevant for future collider facilities, such as FAIR and NICA. By modeling quarkonium at finite velocity we give predictions for the ψ /J/ψ ratio.

The Gauss law potential model
In order to fully utilize the advances in lattice QCD computations of the in-medium heavy quark potential in phenomenological studies, an analytic parametrization of the complex V (r) is required. Such a parametrization may also provide a starting point for modeling heavy quark interactions in regions of the QCD phase diagram currently inaccessible to first principles Monte-Carlo simulations.
More specifically, we require easily evaluable expressions for both the real and imaginary parts of the in-medium heavy-quark potential that provide a faithful reproduction of non-perturbative lattice data where available. In particular it must be applicable in the temperature regime close to and around the crossover transition, where HTL perturbation theory by itself is not valid. In this study we deploy a similar physical reasoning as in [28,30] to construct such an analytic parametrization of the potential, overcoming the previous shortcoming in that we avoid any ad-hoc assumptions on the linear part of the in-medium potential.
The starting point is the fact that the vacuum behavior of quarkonium bound states is described well by the Cornell potential of Eq. (1.4). We then consider this heavy-quark two-body system immersed in a weakly coupled medium described by the HTL permittivity. Both the Coulombic and the string-like part of the Cornell potential will then receive inmedium corrections, which we compute using linear response theory.

Constructing the in-medium model
In linear response theory the electric field at finite temperature, or equivalently the electric potential, can be obtained from its vacuum counterpart by multiplying it by the inverse of the static dielectric constant in momentum space [35], The permittivity, defined as an appropriate limit of the real-time in-medium gluon propagator, imprints the medium effects onto the potential here. In the following we will consider separately the Coulomb and string-like parts of the Cornell potential in the above relation. Note that Eq. (2.1) does not rely on a weak-coupling assumption and remains valid as long as the vacuum field is weak enough to justify the linear approximation. Using the convolution theorem it can be recast in coordinate space as follows, where ' * ' represents the convolution.
To continue we consider the other central building block of our approach, the generalized Gauss law, which holds for electric fields of the form E vac (r) = −∇V vac (r) = qr a−1r . This reduces to the well-known form for Coulombic potentials with a = −1, q =α s , while the linearly rising string case corresponds to a = 1, q = σ. For a general a we have − 1 r a+1 ∇ 2 V vac (r) + 1 + a r a+2 ∇V vac (r) = 4πqδ(r) . (2.4) Denoting the differential operator on the left-hand-side above as G a , we now apply it to Eq. (2.2) to deduce the general integral expressions for each term in the in-medium heavyquark potential: where we have used Eq. (2.4) and that the convolution commutes with G a . For the Coulombic and string cases, respectively, this gives At this point we introduce the explicit expression for the coordinate space in-medium permittivity obtained from the perturbative HTL expression in momentum-space [35], The real part of the real space expression can be calculated via a contour integral and the residue theorem to give while the imaginary part is expressed by a Meijer-G function, (2.10) Let us use Eq. (2.9) and Eq. (2.10) to solve for the in-medium modified Coulombic part of the potential. We find that our ansatz, as expected, reproduces the well known HTL result [16,36] given in Eq. (1.2): The next step is to turn to the string part, Eq. (2.7), for which the formal solution can be straightforwardly written down as The constants c 0 and c 1 will be chosen to ensure the physically-motivated boundary conditions ReV S (r)| r=0 = 0, ImV S (r)| r=0 = 0 and ∂ r ImV S (r)| r=0 = 0. This leads to the following analytical form: 14) With both the Cornell and string in-medium solution at hand we combine the expressions to form the full Gauss law model Let us inspect the properties of each contribution to the in-medium potential found so far. For the real parts, the short distance r → 0 limit, corresponding to each heavy-quark not being able to "see" the intermediate medium, recovers the vacuum Cornell potential. The zero temperature limit corresponds to m D → 0, which also recovers the Cornell potential. At large distances the real part exhibits an exponential flattening-off ∼ e −m D r , which is the characteristic and well-known Debye screening behavior. Due to the extra factor of m D in the denominator, the string contribution to the in-medium real-part will become more and more suppressed at high temperature, eventually giving way to the pure HTL result. The imaginary part arising from the Coulombic contribution asymptotes to a constant at large distances, which is expected for Landau damping. Only the imaginary string part, Eq. (2.15), at first sight appears problematic, as it diverges logarithmically at large r. As we argue in the next section this is a manifestation of the absence of string breaking in the Cornell potential and we will account for it by introducing a well motivated regularization.

Consistent treatment of string breaking
In the preceding section we found that only the imaginary string part, Eq. (2.15), shows an unphysical behavior, in that it diverges logarithmically at large r. To understand the origin of this artefact let us consider the ingredients used. The generalized Gauss law (Eq. (2.3)) is formally correct, and the linear-response relation in Eq. (2.1) is valid under a weak-field ansatz. However, the vacuum potential and in-medium permittivity both operate under assumptions that can be challenged. Firstly, we have utilized a vacuum potential that contains an unending and unphysical linear rise. Secondly, the expression for the complex permittivity given in Eq. (2.8) is a hard-thermal-loop result and as such is strictly only valid at temperatures much larger than T c .
In our computation, both issues manifest themselves in Eq. (2.13), which can be written, after substituting the imaginary part of Eq. (2.8) into Eq. (2.13) and performing the angular integration of the inverse Fourier transform, as follows: (2.17) We have arranged the momentum factors as above to make clear their different origins: the first term (p 2 ) arises from integrating in spherical coordinates and the second (sinc(pr )) after completing the polar integration. The last two terms encapsulate the contribution from the gluon propagator, and it is the 1/p factor here that we identify as causing the weak infrared divergence. In order to regularize this integral, we modify the last term as follows: 1 where ∆ will be a suitably chosen regularization scale. In Eq. (2.17) the spatial integrals can be carried out analytically, which combined with the regularization above gives our new definition of the string imaginary part, where we have also chosen the constant terms to impose the boundary conditions as before.
Eq. (2.19) can be numerically evaluated very quickly, and the only remaining step is now to choose the regularization scale ∆.
To this end we propose the following physically motivated scheme. Note that if we rescale momentum via p → p/m D and rearrange slightly, Eq. (2.19) takes on a suggestive form: where χ(x) = 2 ∞ 0 dp 2 − 2 cos(px) − px sin(px) and ∆ D = ∆/m D . That is, we can express the string imaginary part as some temperature dependent prefactor with dimensions of energy, multiplied by a dimensionless momentum integral. This is identical to the Coulombic expression, where the integral asymptotes to unity in the limit r → ∞. We thus impose the same condition for the string part. This procedure also recovers the correct behavior at large T (large m D ), i.e. the string contribution to the imaginary part diminishes while the Coulombic part grows in stature and we eventually recover the pure HTL result. The value of the regularization parameter ∆ D can be computed numerically (the trigonometric terms in Eq. (2.21) drop out in the r → ∞ limit). Furthermore, since it is expressed in terms of the Debye mass it remains constant and the computation need only be performed once. We find In order to check whether our particular choice of regularization influences the end result we also implemented instead a factor of tanh(p/∆ ) in the integral in Eq. (2.17). The hyperbolic tangent rises linearly at small p and converges exponentially quickly to unity and thus is able to fix the infrared divergence while leaving the ultraviolet behavior unchanged. We find that using this alternative leads to an equivalent ∆ D that lies within 1% of the value in Eq. (2.22), and the subsequent results for ImV are indistinguishable by eye to those obtained via the original method. That the regularisation process is independent of the exact technique used serves as an encouraging cross-check. Fig. 1 depicts the values of the real (left column) and imaginary (right column) part of our novel Gauss law parametrization at three different realistic combinations of temperature and Debye mass. The red solid lines correspond to the Coulombic contributions and the green lines to the string parts. The total is shown as a black solid line. In the panels for ReV the vacuum Cornell potential is added as a gray solid line. For completeness the unregularized string imaginary part is included as a blue line. The figures clearly show how the perturbative HTL results, i.e. the purely Coulombic contribution, dominates at high temperatures.

Comparison to other models
Let us compare the expressions derived above from the Gauss law with other models previously deployed in the study of heavy quarkonium. Before the realization that the in-medium potential is a complex quantity, the focus lay on modeling a real-valued potential. A classic work in this regard is the study by Karsch, Mehr and Satz, who argued based on the twodimensional Schwinger model that the in-medium potential should show both the standard Debye screening term for the Coulombic term as well as an exponential damping factor for the string part: (f) T = 1 GeV, m D 2.25 GeV Figure 1: Temperature dependence of the Coulombic (red), string regularized (green) and string unregularized (blue) imaginary part. In the high temperature limit we recover the purely Coulombic HTL result.
Eq. (2.23) reduces (by construction) to the Cornell potential in the limit m D → 0. At first sight the exponential damping of the string part may appear similar to our result for V S (r). A quantitative inspection however shows that the KMS potential exhibits a stronger dependence on m D , i.e. for the same value of m D the deviation from the m D = 0 limit is more pronounced than in our case.
The KMS potential may also be obtained, as presented in version 1 of [37], by modeling the inter-quark interactions as an effective one-dimensional string like interaction. In addition, the authors postulate that entropy may contribute to the inter-quark interaction. Their heuristic arguments, resting on an identification of the real part of the potential with thermodynamic quantities, leads them to propose for the in-medium string part This expression turns out to be the same as our result, which in this paper has been systematically derived from the Gauss law ansatz. We are thus able to offer an explanation for the presence of the additional r dependent term that does not rely on the ad-hoc identification of the real part of the potential with either the free or internal energies of the quarkonium system. In the context of purely real-valued in-medium potential models, the generalized Gauss law was used for the first time in [31]. In this study a first attempt was made to combine the Gauss law and Debye-Hückel theory to implement the screening of both the Coulombic and string-like parts of the potential. This represented an important step forward towards a systematic modeling of both terms at finite temperature. The authors encountered difficulty in using this parametrization to capture the behavior of the color singlet free energy in lattice QCD, which was taken as a proxy for the in-medium real-part of the potential. This led to the introduction of an additional parameter κ to compensate for these deviations: (2.25) Unfortunately at that time it was not possible to relate κ to either the parameters of the Cornell potential or to the Debye mass. An important step towards parametrizing the inter-quark potential as a complex quantity was taken in [32], where the authors proposed taking the linear response relation in Eq. (2.1) at face value and introduced the same HTL permittivity as used in this paper to directly carry out the inverse Fourier transform. Their computation led to the following proposal for the in-medium potential: where φ was defined in Eq. (1.3) and χ 0 is given by In this approach the modification of the real and imaginary parts are governed by a single temperature dependent parameter, just as in our study. If we inspect the functional form of this parametrization closer, two properties become apparent which challenge the validity of the result. One can be remedied but the other hints at a more foundational difficulty of the approach. Let us first consider the imaginary part. Just as in our derivation, the imaginary string part of the potential diverges logarithmically due to the 1/z term. This divergence may be avoided, as discussed in the previous chapter, by regularizing the unphysical linear rise to infinity of the Cornell potential, i.e. by introducing string breaking. On the other hand if we take a look at the in-medium real part arising from the vacuum string-like potential, we find that it contains an unscreened 1/r term. This would suggest that the deconfined color charges are unable to screen the interactions and a long-range component remains in ReV . Such a behavior is counter-intuitive and does not agree with current lattice data determinations of the in-medium potential.
The authors of [37] recently proposed a very different derivation of the real part in Eq. (2.24). They base it on the direct Fourier transform of the gluon propagator, which however receives an additional non-perturbative contribution originally suggested in [38]. In that construction one also encounters a divergent imaginary part, which may be regularized, as shown in the previous section. The additional term in the gluon propagator is related to a non-vanishing gluon condensate, which has previously been used to justify a a similar real part deployed in the T-matrix approach in [9]. In addition that study models quarkonium screening with different screening masses for the Coulombic m D and string like part m D and a third parameter c s to take into account string breaking effects In [30] the Gauss law was used for the first time in the context of a complex valued in-medium potential. The study brought together the ideas of the Debye-Hückel theory from [31] with the HTL permittivity as used in [32]. However, as already mentioned in the introductory section, to solve the Gauss law equations in that paper, the authors introduced ad-hoc assumptions about the string part of the in-medium potential. The present paper, while using a very similar combination of HTL permittivity and linear response theory, provides a rigorous derivation of the in-medium expressions for ReV (Eq. (2.11) and Eq. (2.14)) and ImV (Eq. (2.12) and Eq. (2.19)) that form the central result of this chapter.

Vetting with lattice QCD data
The most important benchmark for any parametrization of the in-medium heavy quark potential is whether it is able to reproduce the non-perturbative lattice QCD results. Since our Gauss law approach uses the HTL permittivity to modify the non-perturbative vacuum potential it is by no means obvious that it can capture the physics of the inter-quark potential in the non-perturbative regime around the crossover transition. We will show in this section that it indeed works excellently even in this regime.
One hint at why the Gauss law may work where HTL alone is no longer valid, is given by the form of the Gauss law equation of the Coulombic part in coordinate space. As was discussed in [30], using the HTL permittivity leads to an expression that has the same form as a linear response equation for V C as the original Debye Hückel theory, and the Debye mass parameter governs the strength of the linear response. This allows to smoothly connect the expression at finite T to the unscreened potential at T = 0.
The vetting is carried out using published state-of-the-art lattice QCD data for the real-part of the potential [28,29]. The values for ReV and ImV have been extracted from simulations of the HotQCD collaboration on 48 3 × 12 lattices featuring N f = 2 + 1 flavors of dynamical light quarks discretized with the asqtad action. The pion mass on these lattices is larger than physical at m π ≈ 300 MeV, with a transition temperature T c ≈ 175 MeV. As temperature is changed on these lattices by changing the lattice spacing, there are zero temperature ensembles available for calibration purposes.
The first step in applying the Gauss law parametrization consists of fixing the vacuum parametersα s , σ, and c appearing in the Cornell potential, which characterize the test charges inserted into the medium. Two sets of low temperature results for ReV are available [28], to which Eq. (1.4) is fitted. The results are given in Table 1 and shown in gray in Fig. 2. The naive Cornell ansatz works excellently in describing the lattice data in the phenomenologically relevant range of distances 0.1 fm < r < 1.2 fm.
Once the vacuum parameters are determined, the Gauss law parametrization contains Table 1: Best fit result for the vacuum potential parameters after fitting to the low temperature lattice QCD data.
1.781 ± 0.059 2.648 ± 0.042 a single temperature dependent parameter, the Debye mass m D , which controls both ImV and ReV . Here we will fit the values of m D using the real-part. As can be seen from the left panel of Fig. 2 the fit works excellently and the Gauss law parametrization is able to capture the behavior of the potential from the Coulombic region, through the intermediate regime and up to the screening regime at large distances and high temperatures. This gives us confidence that the derived expressions for the in-medium potential do capture the relevant physics encoded nonperturbatively in the lattice data. The best fit values for the Debye mass parameter are listed in Table 2.
Once m D is fixed we can compare the corresponding Gauss law prediction for the imaginary part with the tentative values extracted on the lattice. As shown in the left panel of Fig. 2 we find that the agreement is also excellent down to T = 164 MeV, which on these lattices corresponds already to the confined phase. Only at T = 148 MeV, deep in the hadronic regime, deviations start to appear. Note that the imaginary part in the Gauss law rises more steeply as temperature increases but the asymptotic large distance value behaves non-monotonously with temperature, reflecting the competition between ImV arising from the string and Coulombic parts of the potential.
We have checked that the properly derived Gauss law expression for ReV presented in this paper leads to a overall better fit of the lattice data than the functional form used previously in [28]. We have also checked that an unscreened Coulombic component as proposed in [32] will lead to a worse fit, disfavoring that model.
We conclude that the Gauss law parametrization provides a viable option to summarize the in-medium behavior of the non-perturbative in-medium heavy quark potential using three vacuum parameters and one temperature dependent Debye mass m D .

Extension to a running coupling
In anticipation of upcoming high resolution lattice QCD computations of the in-medium heavy quark potential, it is prudent to consider the effects of a running coupling in the Gauss law parametrization. While in the simulation data deployed in the previous section the short distance regime was still well described by a naive Cornell potential, more recent lattice studies of heavy quark interactions [39] have shown that at shorter resolved distances the running will manifest itself. Thus we consider the strong coupling parameter of our Cornell potential to become a function of distanceα s →α s (r) and writẽ Note that in the context of the vacuum potential in Eq. (1.4), we have already implicitly included the termsα (1) s andα (2) s by absorbing them into the other vacuum parameters: In a thermal setting, this would necessitate including r a terms other than a = −1, 1 in the formulation of the in-medium potential.
To do this, we must use the generalized Gauss law operator G a given in the left-handside of Eq. (2.4), but with a modified right-hand-side that includes the real-space complex permittivity (following the procedure in Sec. 2.1) (2.32) With the real space expressions given in Eqs. (2.9) and (2.10), a computer algebra program such as Mathematica will give a general solution for general a as follows: where G is again the Meijer-G function and Γ is the upper incomplete Gamma function defined via Note that the imaginary part of the HTL solution expressed via the integral Eq. (1.3) is equivalent to Eq. (2.34) with a = 1, as noted by the authors in [16]. Two remarks are in order: firstly, there is nothing in principle that prohibits the expressions above from being used throughout the remaining analysis in this chapter. We have decided not to do so because currently published lattice data on the in-medium potential does not yet reach the regime where the running is significant. Secondly, for a ≤ 1 the imaginary part exhibits a divergence. In order to employ Eq. (2.34) in a phenomenological study one would need to carry out a regularization procedure similar to that discussed in the last section, which is in principle possible but we have not investigated this further.
Keeping in mind that the Gauss law approach can straightforwardly be extended to accommodate a running coupling we nevertheless proceed to use its naive formulation (a = −1, 1 only) in order to investigate the in-medium properties of heavy quarkonium in subsequent chapters.

Extension to finite velocity
In preparation for the study of heavy quarkonium at finite transverse momentum we need to consider how to extend the Gauss law model to treat a heavy quarkonium bound state moving through the plasma at finite velocity. To this end we will follow the ideas laid out in [40]. That is, one considers a QCD plasma in thermal equilibrium and a reference frame in which the medium moves at velocity v with respect to the bound state at rest, a so called hot wind scenario. It then becomes necessary to distinguish two separate alignments: one in which the medium velocity is parallel to the dipole axis of the bound state and another in which it is perpendicular. This leads to distinct self-energies with different angular dependencies and correspondingly, different expressions for the potential in each of the alignments. As reviewed in detail in Appendix A the corresponding in-medium permittivity can be computed for both alignments and used to set up a finite-velocity Gauss law model in the same manner as in Sec. 2.1.
For the parallel alignment case, we find for the Coulombic part (2.36) and The retarded self-energy Π R is given in Eq.
where Π ⊥ R is given in Eq. (A.15). For the string part, we invite the reader to consult Appendix A. The method is identical to the static case in the sense of the integration in Eq. (2.13), now with a modified finitevelocity permittivity. Since the resulting expressions offer no further intuition we omit them here, but note that the constants are again set to impose the same physical boundary conditions as in the static case.
In Fig. 3 we show a selection of finite-velocity potentials at T = 155 MeV both in the parallel and perpendicular alignments. Let us focus on the real part first. The velocity dependence for the two alignments is significantly different. While for the case r v the real part is weakened with increasing v, the opposite happens for r ⊥ v and one eventually recovers the Cornell potential. As had been pointed out in [41] the reason for this non-trivial behavior is the fact that the quarkonium state encounters a different effective temperature depending on the magnitude and orientation of the velocity relative to the medium. The higher the velocity, the more the effective temperature deviates from its v = 0 value, with a hotter region in the forward direction and cooler region in the backward direction. Since in the permittivity one integrates over the orientation θ, it seems that for r v the hotter region dominates, while for r ⊥ v the colder region contributes more significantly. The fact that the real part becomes negative at large distances is not troubling, as it simply tells us that similar to the centrifugal barrier present in P-wave states, there is now a finite probability of the in-medium S-wave tunneling into an unbound configuration. The r v imaginary part of the potential behaves in a well defined manner up to v = 0.9c, in that it shows non-monotonicity but stays positive at all distances. Once ImV becomes negative (in the sign convention used here) it may in principle introduce an instability when being used in a Schrödinger equation, rendering the computation unreliable. Even though the potential does show excursions of that manner for v > 0.9c, in practice we have not encountered any numerical difficulty when solving for the corresponding spectra presented in Sec. 4.3. For r ⊥ v the imaginary part diminishes at higher velocity, as if the quarkonium state encounters a cooler and cooler surrounding, consistent with the real part in that scenario.
In the subsequent chapters we will explore how this non-trivial modification of the potential translates into changes in the in-medium spectral functions and attempt to gain a first glimpse into production yields of heavy quarkonium in heavy-ion collisions at finite transverse momentum.
We would like to note that in order to comprehensively capture the physics of heavy quarkonium with a finite center of momentum motion, a genuine non-equilibrium realtime approach eventually has to be developed. Efforts in this direction in the context of the open-quantum-system approach to heavy quarkonium are ongoing and will eventually incorporate information from the Gauss law model developed here.
3 Quarkonium spectral properties at finite temperature

Continuum corrections
Before we can embark upon studying the in-medium properties of heavy quarkonium based on the lattice vetted Gauss law model, we have to acknowledge the fact that the lattice data used in the previous section was not continuum extrapolated. While there is activity in the community to bring the extraction of the potential closer to the continuum limit [42], no truly extrapolated results are available today. Thus we need to manually correct for the discrepancy between the discrete lattice results and those that eventually will lead to agreement with experiment.
To this end we follow the strategy laid out in [28]. The idea is first to determine a phenomenological set of vacuum parameters that, via solving the Schrödinger equation, reproduce the masses of the known quarkonium ground state particles. In addition we will then use an appropriately rescaled version of the lattice-fitted Debye mass to implement the finite temperature effects. The starting point is the bottomonium system, which due to the large mass of the bottom quark is the most amenable to the potential description. Furthermore, since the bottom mass is much larger than Λ QCD , the matching of pNRQCD with QCD can be carried out perturbatively in vacuum and one finds that the appropriate mass to use in the Schrödinger equation is the so called renormalon subtracted mass [43] m RS b = 4.882 ± 0.041 GeV. are able to reproduce both the S-wave and P-states very well (see Table 3). The next step is to consistently determine the only remaining unknown parameter-the charm quark mass. Since the heavy quark potential is a universal quantity in the sense that at lowest order in pNRQCD the same expression is used for both heavy quark families, we expect the vacuum values in Eq. (3.2) to remain the same for charmonium. Thus by reverse engineering the procedure used above, we can 'fit' the charm mass to reproduce the lowest S-wave states (J/ψ, ψ ) for our vacuum parameters. The best-fit value reads m fit c = 1.4692 GeV. Since finite mass corrections are more important for charmonium, we expect the agreement to be worse in this case, which is indeed observed in Table 4. All necessary parameters required for describing the T = 0 physics of quarkonium within the potential approach are now set. The next step is to consider how the lattice discretization affects the fitted values of the Debye mass parameter m D . In lattice calculations the light quark masses do not take on their physical values and the chiral crossover temperature is increased. For the lattice spacings used in this study, T lat c = 172.5 GeV. The goal is now to undo the difference between this and the physical crossover temperature by an appropriate rescaling. To this end we consider a dimensionless ratio of the lattice fitted m D (t) / σ(t) where t = T /T c . The slight dependence of the lattice string tension on the lattice spacing is taken into account in this ratio. The continuum corrected value of the Debye mass is then taken to be m phys where m D (t)/ σ(t) is given in Table 1 and σ cont from Eq. (3.2). In order to explore the changes in the heavy quark potential at different temperatures, it is useful to also parametrize the temperature dependence of the Debye mass itself. Let us first look at the perturbative expression for the Debye mass. With dynamical quark masses m u,d set to zero, the leading order result for SU (N c ) with N f fermions at zero baryon chemical potential is [45] m phys where Λ = 2πT is the renormalization scale. It is well established that the Debye mass can only be calculated up to leading order plus a logarithmic correction at next to leading order before truly non-perturbative contributions come into play [46]. For the purposes of this study, we account for this via two additional terms m phys The non-perturbative constants κ 1 , κ 2 will be fixed by performing a fit against the continuum corrected lattice results Eq.  The resulting plot of Eq. (3.6) is shown in Fig. 4 together with the continuum corrected Debye mass points. Note that the perturbative leading order result in Eq. (3.5) leads to a increase in m D /T when we approach T c from above. As we eventually need to recover the T = 0 Cornell potential below T c the true behavior needs to exhibit a downward trend eventually, as it does in the lattice determination. This deviation from the perturbative behavior is easily captured by κ 1 and κ 2 .

In-medium spectral functions
Spectral functions provide the quantum field theoretical answers to questions about particle properties, i.e. their masses and decay widths. Our goal is to learn about the in-medium properties of quarkonium from an inspection of the thermal spectral functions computed from the Gauss law potential. To this end we follow the Fourier space method introduced in [48], where the following Schrödinger equation is established in describing the time evolution of the vector channel unequal-time point-split meson-meson correlator D > (t; r, r ):  A detailed discussion on how to explicitly compute the spectral functions from Eq. (3.8), both formally and practically, can be found in Appendix A of [48]. Furthermore, the pseudoscalar and axial vector channels were shown not to include any qualitatively new structures, with In this study we focus on the S-wave quarkonium spectral functions and show in Fig functions, which exhibit three well defined bound states for bottomonium and two for charmonium below threshold. Qualitatively similar to what has been reported in the literature [28], we find characteristic changes as temperature increases. Both a broadening of the peaks, as well as a shift of their central value to smaller frequencies is observed. We note that the strength of the change is clearly ordered with the vacuum binding energy of each individual state, where E bind is defined from the distance between the spectral peak and the threshold. Just as intuition predicts, the more deeply the state is bound the less susceptible it is to medium effects. This behavior is observed consistently in potential based computations and is also in agreement with recent studies of directly reconstructed bottomonium and charmonium spectral functions from lattice QCD [14].
To more quantitatively explore the in-medium properties we can consult scattering theory. If a narrow resonance pole lies close to the real frequency axis then its spectrum can be described by a Breit-Wigner distribution. On the other hand, if it features a significant decay width we may employ a skewed Breit-Wigner of the form which is able to disentangle the bound state signal from the background continuum. Here, E denotes the energy of the resonance, Γ its width and δ the phase shift. The constant terms C 1 and C 2 model artefacts beyond the spectral peak we are interested in.
For completeness we plot in Fig. 6 and Fig. 7 the in-medium masses, as well as the threshold behavior for charmonium and bottomonium respectively. The results are consistent with previously obtained quarkonium properties based on the previous formulation of the Gauss law model in [30].

Applications to Heavy Ion Collisions
The computation of the in-medium spectral functions has already provided us with vital insight on the properties of heavy quarkonium in thermal equilibrium. The question however remains of how to connect this information to actual measurements carried out in heavyion collision experiments. Contrary to light mesons, where a direct link exists between in-medium spectral functions and the measured decay leptons, for heavy quarkonium we do not measure their thermal decay. Instead, at hadronization a number of vacuum states are created whose decay is recorded long after the QGP has ceased to exist. Thus we need to translate all in-medium information into abundances of vacuum states to connect to experiment.
In previous studies the Gauss law model has been deployed to investigate the production of quarkonium in heavy-ion collisions in the simplest possible scenario, where the quarks are almost at rest with the surrounding medium. In addition, since lattice QCD simulations have not yet produced results for the heavy-quark potential at finite baryon density, only predictions for highest energy LHC collisions have been presented.
Here we will go beyond these results by modeling heavy quarkonium production both at finite baryo-chemical potential, relevant for lower-energy heavy-ion collisions at RHIC, FAIR, and NICA, and at finite transverse momentum.

ψ to J/ψ production ratio
Among the currently most highly sought observables at the LHC heavy-ion program is the ratio of ψ to J/ψ produced in Pb-Pb collisions. In contrast to the nuclear modification factor R AA for each individual species, this ratio promises to be highly discriminatory between different phenomenological models for quarkonium production in heavy-ion collisions. Our goal is therefore to estimate this ratio.
Following the ideas laid out in [28], we utilize the in-medium spectral functions computed in the preceding sections and will assume kinetic thermalization of the charm quarks in the late stages of a collision.. The idea is to convert in a meaningful fashion the in-medium spectral information about ψ and J/ψ into the number of produced vacuum states. The assumption behind this step is that of an instantaneous freeze out, where at T = T c the in-medium particles abruptly change into vacuum particles.
To translate in-medium spectral peaks into vacuum states we consider the relation between the dilepton emission rate and the in-medium spectral function [49]: (4.1) Here, n B denotes the Bose-Einstein factor, Q q the electric charge of the heavy quark in units of e, and α e the QED coupling constant. The four-momentum is denoted P = (p 0 , p) and the finite mass of the leptons has been neglected. Eq. (4.1) relates the weighted area under the in-medium peak to the rate of dileptons emitted from that state. Thus, integrating the right-hand-side above gives the number of in-medium lepton pairs produced by each of the states. The prefactors drop out if only the ratio is computed, leaving Let us emphasize again that this quantity is not what is measured in experiment. Since the plasma is diluted away long before the charmonium states decay, the number of dileptons eventually measured originate from the vacuum-state remnants of the in-medium structures observed here.
Thus we must project the states corresponding to the finite temperature peaks onto the T = 0 vacuum states. The computation proceeds as follows. At leading order the vector channel spectral function ρ V depends only on P 2 = p 2 0 − p 2 ; after performing a change of variables to ω = p 2 0 − p 2 , Eq. (4.2) becomes In this expression, the contribution from each bound state arises from the corresponding peak area in ρ V (ω) /ω 2 . We fit each peak structure with the skewed Breit-Wigner in Eq. (3.13), thus allowing the different contributions to be distinguished and the thermal mass M n and width of each to be ascertained. Now, the projection onto the vacuum states is implemented by writing ρ V (ω) /ω 2 = A n δ(ω − M n ); that is, we allow the in-medium states to collapse into delta peaks that represent the vacuum states, while retaining the peak area A n to account for the different contributions. We have also confirmed numerically that it is indeed possible to approximate the Breit-Wigner peak by a delta function in this manner. Imposing this on Eq. (4.3) and carrying out the now trivial ω integral then gives Switching to spherical coordinates in momentum then leaves only a dp integral that can be easily evaluated numerically. Finally, in order to obtain the total number density we must divide by the electromagnetic decay rate of the vacuum state, which is proportional to the square of the wavefunction at r = 0 divided by the square of the mass of the state [50]. These values we calculate (4.5)

Finite µ B phenomenology
In order to make use of Eq. (4.5) in connecting to experiment, we require a prescription to evaluate our Gauss law potential model (and hence the resulting spectral functions) at a given centre-of-mass energy. The strategy here is two-fold. Firstly, we note that the successful statistical hadronization model provides a well-established scheme with which to estimate the thermal parameters of the produced medium at chemical freeze-out resulting from a collision at a given √ s N N . Starting from the grand canonical partition function of known hadrons, one is able to reproduce the measured yields by adjusting the three parameters of the model: collision volume, temperature, and baryo-chemical potential (µ B ). The most recent values [51] are: (4.6) where √ s N N is the dimensionless numerical value of the centre-of-mass energy measured in GeV. These are plotted in Fig. 8. In this model the freeze-out temperature quickly asymptotes to the limiting temperature of 158 MeV while µ B drops monotonously to almost zero at high collision energies such as those probed at RHIC and LHC. NA50 Pb-Pb, √ s N N = 0.017 TeV, central, 0 < y < 1, p T > 0 GeV ALICE Pb-Pb, √ s N N = 0.017 TeV, 0-20%, 2.5 < y < 4, 0 < p T < 8 GeV CMS Pb-Pb, √ s N N = 2.76 TeV, 0-100%, |y| < 1.6, 6.5 < p T < 30 GeV CMS Pb-Pb, √ s N N = 5.02 TeV, 0-100%, |y| < 1.6, 6.5 < p T < 30 GeV Figure 9: The prediction of this work (green) for the relative production yield of ψ to J/ψ. We also include the statistical hadronisation model prediction [51] (purple) and experimental data measured by the NA50 [52], ALICE [53] and CMS [54,55] collaborations (red) for Pb-Pb collisions as well as the pp baseline [51,56] (orange).
Secondly, note that the medium effects in our potential model are captured entirely by the value of the Debye mass m D . Eq. (3.4) gives our interpolated continuum-corrected expression at µ B = 0, denoted now as m D (T, µ B = 0). We now postulate how to extend that formula into the realm of finite baryo-chemical potential. At leading order, the Debye mass can be directly calculated perturbatively at finite baryo-chemical potential [45]. We propose to add this µ B -dependence to the temperature-only dependence of Eq. (3.4): Here, the renormalization scale is also modified to Λ = 2π T 2 + µ 2 B /π 2 . This expression is valid only for small values of µ B . On the other hand, at very large values of µ B the chemical potential itself becomes the only relevant scale and we expect from dimensional grounds that such a scale enters in m D again linearly. Thus our modeling assumption is to naively adopt Eq. (4.7) at all values of µ B to extend the Gauss law parametrisation into the finite baryo-chemical potential regime. In the absence of reliable lattice data at non-zero quark chemical potential, we hold the non-perturbative constants κ 1 and κ 2 in m D (T, 0) the same and include the chemical potential dependence solely by adding the linear part as shown above. Under these assumptions we may explore the domain of finite µ B and via the information provided from the statistical model, also at different √ s N N (i.e. via Eq. (4.6)).
With all ingredients in place, we may now proceed with calculating the ψ to J/ψ ratio over a range of center-of-mass energies. The results from this entire procedure are plotted in Fig. 9, as well as a comparison with the statistical hadronization model prediction and the most up-to-date experimental data for Pb-Pb collisions and the pp baseline. This extends previously available computations, valid only at the highest beam energies, to those relevant for NICA and FAIR.
Our analysis is based on in-medium spectral functions and is independent from that performed by the statistical hadronization model. The only information shared among the two are the values for T and µ B extracted from the yields of light hadrons. We find as expected from the previous Gauss law studies that the results at vanishing µ B lie very close to the prediction from the statistical model. In the lattice fits carried out using the new and improved Gauss law model we have been more conservative in the estimation of our uncertainties, which is why the present results are fully compatible with the statistical model. Note that while different in form, both approaches share that they consider a fully kinetically thermalized scenario. A good agreement between the two and the experimental data thus further supports the interpretation that charmonium at LHC has reached a significant degree of kinetic equilibration with its surrounding.
We find that extending our Gauss law model to finite µ B , i.e. lower √ s N N , the agreement with the statistical model persists. Even though our assumptions to do so were rather crude, they lead to both a qualitatively and even quantitatively very similar trend for ψ /J/ψ. On the other hand the full validity of these results when compared to experimentally measured data is somewhat questionable. Whether charmonium can be considered as kinetically thermalized in collisions as low as √ s N N ∼ 40 MeV, for example as carried out at RHIC, remains to be seen, in particular in light of the difficulties of measuring a finite elliptic flow for J/ψ there.

Finite transverse momentum
Not long after the initial interest in quarkonium as a probe of the QGP, the first works began to appear that considered how the naive Coulombic Debye screening description could be extended to account for quarkonium moving with a finite velocity [57]. Recent years have seen a revival of interest in this direction, with new approaches employing modern effective field theory techniques to tackle the problem [40,41,58,59]. The phenomenological motivations are clear; quarkonia produced in heavy-ion collisions traverse the hot medium before being measured with finite transverse momenta p T , and accounting for this may lead to qualitatively new QGP phenomena such as the formation of wakes [60][61][62].
In order to take a first step towards realistic phenomenology based on the finite velocity Gauss law model constructed in Sec. 2.6, we require a prescription for converting from a general medium velocity v to the transverse momenta p T commonly measured in heavy ion collisions. Such a prescription has been described in [59], which we briefly review here. Consider a heavy quarkonium traversing the QGP with a momentum P µ = p 0 , p measured in the lab frame. The QGP will also have a velocity in that frame, denoted w, and it is the relative motion that needs to be estimated and subsequently employed as v in the Gauss law potential. Assuming a central collision and constant w throughout, a typical value for the LHC is w ∼ w ⊥ ∼ 0.66. Further assuming a thermalized and isotropic system, the modulus of the relative velocity will be given as where M is the heavy quarkonium mass. Note that this depends on the angle ϕ between the QGP velocity w and the quarkonium 3-momentum p.
Following the same procedure as the preceding sections, we may now calculate finite velocity spectral functions. Some representative examples are shown for charmonium in Fig. 10. The qualitative difference in the potential for parallel (top) and perpendicular (bottom) alignments manifests itself also in the the behavior of the spectral function. In the perpendicular case, going to higher velocities is reminiscent of increasing temperatures in that the peak, generally speaking, is broadened and shifted to lower frequencies. Note however that this trend is eventually reversed at ultra-relativistic velocities. When investigating the physics of bottomonium at finite velocity in an EFT picture [41], a similar phenomenon was encountered with increasing velocity leading to spectral modifications similar to an increased temperature. In contrast, the parallel case exhibits a shift to higher frequencies with little or no effect on the width.
We may also use the procedure of the preceding sections to estimate the ψ to /J/ψ production ratio at finite velocity and through Eq. (4.8) can relate this to a simplified description of a heavy ion collision. In Eq. (4.8) we replace p with p T and further assume a uniform distribution over ϕ (see left panel of Fig. 11) before taking the mean. The results of this method at large center-of-mass energies are shown in the right panel of Fig. 11. We find that the effect of the finite center of momentum motion on the quarkonium state in the fully thermal Gauss law model is moderate. Up to the p T = 25 GeV considered here, we find an 11% increase in the production ratio for a parallel alignment of the dipole and up to 17% for a perpendicular alignment. We have to keep in mind that the values obtained here rely on many simplifying assumptions, for example that the expansion of the fireball was taken to be isotropic with a constant velocity. What may help us is that if the production of quarkonium really is dominated by the physics around the freezeout, as suggested by the statistical model of hadronization, then we indeed only need to track the shells of the QCD medium, which at a given moment are close to T = T c . How well the model presented here captures the physics of quarkonium in a heavy ion collision will be testable at the upcoming Run 3 at the LHC, where the first accurate experimental values for ψ /J/ψ are expected to become available.
Note that a recent estimate [63] of the p T dependence of the ratio, combining the statistical model with a more realistic evolution of the temperature profile in the collision center and careful separation of the core and corona of the collision, predicts a strong change up to a factor of 3-4.

Conclusion
Heavy quarkonium is a vital tool in developing our understanding of strongly interacting matter and connecting experimental results to the fundamental theory of QCD. On one hand, progress in effective field theories has demonstrated rigorously the validity of describing their in-medium behavior with an effective potential in an appropriately defined Schrödinger equation; on the other, advancements in lattice techniques continue to provide non-perturbative results against which such approximations can be checked.
The main conceptual result of this work is a rigorously derived model for the in-medium heavy quark potential based on the generalized Gauss law in linear response theory. It describes the in medium modification of the vacuum Cornell potential by self-consistently incorporating a weakly-coupled medium described by the HTL permittivity. The resulting analytic expressions depend only on a single temperature dependent parameter and are able to reproduce the lattice results for the real part of the potential even in the non-perturbative regime close to T c . The string imaginary part showed an unphysical logarithmic divergence, which we attribute to the equally unphysical unending linear rise of the vacuum Cornell potential. By considering the presence of string breaking it is possible to regularize this artefact and we were able to give physically sound predictions for the imaginary part of the potential that qualitatively matched the lattice data. The presented work has improved on the conceptual clarity and technical robustness of the Gauss law model compared to other potential models described in the literature. In preparation of upcoming high resolution lattice QCD data for the inter-quark potential, the straightforward extension of the Gauss law model with a running coupling has been discussed.
Using a combination of weak-coupling computations and dimensional analysis, we introduced an extension of the Gauss law model to finite baryo-chemical potential. The extension to quarkonium moving relative to the QCD medium required us to consider two separate alignments of the quark anti-quark dipole with respect to the velocity when computing the corresponding in-medium permittivity. The resulting expressions for the in-medium potential, while lengthy, could be provided in explicit form. Both extensions of the model are required to step towards describing phenomenologically relevant scenarios in heavy-ion collisions, which are not yet amenable to direct lattice QCD simulations.
A continuum correction on the vacuum parameters, as well as m D , was performed that allowed physically realistic spectral functions to be computed in the S-wave channel, which formed the basis of our phenomenological investigation. Similar to previous studies based on the previous Gauss law model, we find the characteristic broadening and shifting of spectral features to lower frequencies with increasing temperature. The strength of the inmedium modification is hierarchically ordered with the vacuum binding energy. A skewed Breit-Wigner was fitted to each resonance peak in order to obtain quantitative results for in-medium masses and thermal width, which are consistent with previous studies.
The first phenomenological result of this work lies in extending the calculation of the ψ to J/ψ production yield to finite baryo-chemical potential and subsequently to lower beam energies relevant for NICA and FAIR. By assuming an instant freeze-out at around the chiral crossover temperature, we found an excellent agreement with the statistical hadronization model, and our prediction aligned with the latest results from ALICE and CMS to within the experimental errors. As our approach, based on in-medium spectral functions, is largely independent of the statistical model of hadronization but shares the idea of a fully kinetically thermalized quarkonium, the agreement corroborates the interpretation of charm quarks becoming equilibrated in the hot fireball before transitioning into vacuum states at the freeze-out boundary.
The second phenomenological result is our estimate of the change in the ψ to J/ψ production yield for finite transverse momentum. We find increases between 11% − 17% for an increase in p T from zero to 25 GeV, which is moderate compared to predictions based on the statistical model which foresees an increase by a factor 3-4.
We are confident that the availability of this cleanly derived and lattice-vetted complex potential model will be of use to the quarkonium phenomenology community. The Gauss law model described here is future proof, as it is ready to accommodate the upcoming high precision and high resolution lattice data on the inter-quark potential, where for example a running coupling will be relevant. While in this study we were only able to access the information contained in thermal in-medium spectral functions, we are looking forward to seeing the complex potential inform simulations in the open-quantum-systems framework for heavy quarkonium, where a more detailed analysis of the in-medium real-time evolution and recombination dynamics at freeze out are possible.

A Heavy quark potential at finite velocity
The work of [59] first investigated non-relativistic QED bound states moving with finite velocity in a background thermal bath via a rigorous effective field theory approach. This was then generalized to QCD and heavy quarkonium in [41], which forms the starting point of our analysis. The general framework is as follows. We assume that the QCD plasma is in thermal equilibrium at temperature T and consider a reference frame in which the medium moves with velocity v with respect to the heavy quark bound state (QQ) at rest. This frame has been used successfully in the past [64]. The particle distribution functions are given by where the plus (minus) sign refers to bosons (fermions) and Here, γ = 1/ √ 1 − v 2 is the Lorentz factor with v = |v| and the four-momentum is P = (p 0 , p). The study of a bound state in a moving thermal medium is equivalent to studying a bound state in non-equilibrium field theory [65]; in such a formalism the Bose-Einstein or Fermi-Dirac distributions are generalized, which in our case are given by the boosted versions in Eq. (A.1). For a thermal medium of massless (anti-)particles, non-equilibrium field theory gives where p = |p| and θ denotes the angle between p and v. The distribution functions in Eq. (A.1) then become where the effective temperature is defined .
Intuitively, Eq. (A.5) can be understood as a Doppler effect. For v 1 it is shown in [41] that T ef f ∼ T for all directions θ, while for v ∼ 1 the temperature felt by the bound state varies more significantly. The new scales introduced by considering a thermal medium can be understood via light-cone coordinates by defining a maximum (T + ) and minimum (T − ) measurable temperature with T − < T < T + . The subsequent discussion assumes T − ∼ T + , which is strictly not true as v → 1 however in [59] it was found that correct results were obtained for QED by a simple extrapolation.
The authors in [41] then proceed with an inspection of the hierarchy of scales in the formalism outlined above. The details will not be included here, however the argumentation follows the construction of pNRQCD HTL . In the regime T 1/r m D E binding , one may employ the HTL real-time formalism and extend the computation of the heavy quark potential at finite temperature to include the presence of a moving thermal bath. This was first performed for the real Coulombic part in [57,61] before being extended to the newly-understood imaginary part in [59]. 2 More recently, these results were combined with the linear response ansatz in order to model the in-medium and finite-velocity modifications to the string as well as the Coulombic part of the Cornell potential [40].
We will now briefly review the real-time HTL calculation, following the procedure in [40], before applying our model prescription. The physical component of the longitudinal component of the gluon propagator can be written in terms of the corresponding retarded, advanced, and symmetric propagators: where each propagator can be obtained from its corresponding self-energy. Note that we have made explicit the velocity dependence, since these represent different quantities to those discussed previously. In this framework, the symmetric propagator is given as where Π R(S) (p, u) is the retarded (symmetric) self-energy. In the frame where the bound state is at rest, the retarded self-energy can be parametrized as The QQ dipole is always considered to lie along the z-axis, which gives rise to two separate alignments. Firstly, if the velocity direction is parallel to the axis of the dipole, Eq. (A.11) becomes z = v cos(θ) where θ is now simply the polar angle of the momentum vector. One should keep in mind that the momenta here are associated with the mediating gluons and thus align with the dipole direction. In the second case the velocity of the medium lies in the x − y plane and makes an angle β with the x-axis. Labeling the azimuthal angle of the momentum vector as φ, Eq. (A.11) then becomes After some manipulations, Eqs. (A.8)-(A.13) can be combined to give the complex retarded gluon self-energy in the parallel and perpendicular directions: With these expressions, the retarded propagator is obtained via .
(A. 16) Furthermore, the advanced self-energy can be obtained with the relation Similarly, one can calculate the symmetric self-energy for both cases [59]. The result is We now have all of the ingredients to assemble Eq. and for the perpendicular case as D ⊥ S (p, u) = −2πim 2 D T 1 − v 2 3/2 2 + v 2 − v 2 sin 2 (θ) cos 2 (φ − β) 2p 1 − v 2 + v 2 sin 2 (θ) cos 2 (φ − β) 5/2 p 2 + Π ⊥ R (p, u) p 2 + Π ⊥ R (p, u) * . One notices that other than the trigonometric angular factors, the momentum structure takes a very similar form to that given in Eq. (2.8), with the Debye mass being replaced by the retarded self-energies. Indeed it is easily checked that taking the v → 0 limit recovers the static expression.
With the entire framework now in place, we can apply our procedure of modeling the potential via the generalized Gauss law and linear response ansatz. This is again where our analysis takes a different path from the existing literature. As for the static case, our method requires solving two ordinary differential equations (Eqs. (2.6) and (2.7)) but now with a modified right-hand-side that includes the in-medium finite-velocity complex permittivity. The first step is to ascertain the real-space expression ε −1 (r, u). This proves to be somewhat trickier than in the static case, due to the inherited angular dependence of the self-energies. For the parallel case, the integrand for the real part can be manipulated to Re ε −1 (r, u) ∼ π/2 0 dθ ∞ −∞ dp p 2 sin(θ) h(θ, u) Re p 2 p 2 + Π R (p, u) e ipr cos(θ) , (A. 25) where h(θ, u) contain the appropriate angular factors. A similar expression exists for the perpendicular case, with an extra integral over the azimuthal angle. Our attempts to compute this integral via the usual contour techniques led to an ill-defined limit. Thus we instead propose to follow the steps in [40], that is, we solve for the Coulombic part of the potential in momentum space before Fourier transforming back to attain the real space expression. This is formally correct and not conflict with any other part of our procedure. We then assume that the deduced ReV C satisfies our defining in-medium equation, − ∇ 2 V C (r, u) = 4πα s ε −1 (r, u) , (A. 26) and apply the derivatives in order to acquire the in-medium permittivity to be used for the string part. This procedure is further justified since the resulting real part of the real-space permittivity reduces to the correct expression in the v → 0 limit. The final expression for the parallel case, also given in the main text, is For the perpendicular case the results are identical after the replacements in Eq. (2.38). These expressions merit some discussion. Firstly, we highlight that the potentials here differ slightly from those in [40]. In the computation of Eq. (A.27) above, one is faced with a momentum integral like ReV C ∼ d 3 p h(θ, u) Re 1 p 2 + Π R (θ) (A. 29) where the self-energy contains an angular dependence. The strategy is first to exploit a symmetry present in Π R , namely that it is symmetric around θ = π/2. After some manipulations this allows the dp integral to be extended over the entire real domain, such that a contour integration can be performed by continuing p into the complex plane. One can see in the denominator in Eq. (A.29) that the poles exist at ±i √ Π R . Thus taking the real part-in accordance with the retarded propagator definition-must be done after the identification of the pole location. In [40] the taking of the real part was performed on the entire result of the contour integral, which can be seen in the discrepancy between Eqs. (34) and (35) in that text. We have also checked this numerically and confirmed that our expressions are correct. The second rather technical detail that arises from performing the contour integral is that closing in the upper-half plane gives a non-vanishing contribution. This leads to the term ∼ 1/r in Eq. (A.27) above that is necessary to ensure the correct expression is recovered in the limit v → 0.
The computation of the imaginary part of the real-space permittivity is somewhat simpler. It amounts to inverse Fourier transforming Eqs. (A.21) and (A.22) with an extra factor of p 2 . The momentum integrals can be performed analytically and the resulting expressions are angular integrals as we have just seen. In practice we have found that carrying out the same procedure as for the real part, i.e. computing Im ε −1 (r, u) by taking appropriate derivatives of the existing ImV C in [40], leads to an expression that gives the same result as the direct transform, however is numerically more stable and faster to evaluate. The final result that we use for the parallel case is where Shi(x) and Chi(x) are defined respectively by Chi(x) = γ E + log(x) + The perpendicular alignment is again obtained by the replacements in Eq. (2.38). We do not include the corresponding expression for Im ε −1 (r, u) since it is rather long and does not provide any intuition. With a computer algebra program such as Mathematica, it can be easily obtained by acting the Laplacian on Eq. (A.30) or the perpendicular case equivalent. Finally, our expressions for the string part in-medium finite-velocity potential are then achieved in the same manner as for the static case, that is via V S (r, u) = c 0 + c 1 r − 4πσ