---
abstract: 'Based on the observation that the heavy quark - anti-quark potential value at infinity corresponds to twice the D meson mass, we constrain the asymptotic value of the heavy quark potential in a hot medium through a QCD sum rule calculation of the D meson at finite temperature. We find that to correctly reproduce the QCD sum rule results as well as a recent model calculation for the D meson mass near the critical temperature, the heavy quark potential should be composed mostly of the free energy with an addition of a small but non-trivial fraction of the internal energy. Combined with a previous study comparing potential model results for the $J/\psi$ to a QCD sum rule calculation, we conclude that the composition of the effective heavy quark potential should depend on the inter-quark distance. Namely, the potential is dominated by the free energy at short distance, while at larger separation, it has a fraction of about 20% of internal energy.'
author:
- Philipp Gubler
- Taesoo Song
- Su Houng Lee
title: $D$ meson mass and heavy quark potential at finite temperature
---

Introduction
============

Quarkonium is a colorless and flavorless bound state of a heavy quark and its anti-quark. Since Matsui and Satz proposed the suppression of quarkonium in relativistic heavy-ion collisions as a signature of the quark-gluon plasma formation due to color screening [@Matsui:1986dk], this phenomenon has been investigated in numerous theoretical and experimental studies (see, for instance, the recent review in Ref.[@Rothkopf:2019ipj]). One of the most critical quantities affecting quarkonium suppression is the potential between a heavy quark and anti-quark. In lattice quantum chromodynamics (QCD), the free energy of the static quark antiquark pair can be obtained through the rectangular Polyakov loop along the space-time direction. As one increases the temperature, the long distance part of the free energy exhibits a sudden saturation near the critical temperature, which can be interpreted as the onset of deconfinement. If the free energy is adopted for the potential of the heavy quark system [@Satz:2015jsa], this behaviour leads to a weakening of the quarkonium binding with increasing temperature so that the $J/\psi$, 1S bound state of a charm and anti-charm quark, dissolves around 1.1 $T_c$ [@Lee:2013dca]. On the other hand, at finite temperature the heavy quark and anti-quark can themselves be polarized acquiring additional energy. This leads to the internal energy of the system as another candidate for the heavy quark potential, which has an additional contribution from the entropy density [@Kaczmarek:2002mc], $$\begin{aligned}
U=F+TS,\end{aligned}$$ where $F$, $T$ and $S$ are the free energy, temperature and the entropy density $S=-\partial F/\partial T$, respectively. While the free energy at large separation drops rapidly near the critical temperature, the entropy density becomes large, thus acting as a high potential wall, which prevents the quarkonium from dissolving. As a result, if the internal energy is adopted for the heavy quark potential, $J/\psi$ binding persists beyond 1.6 $T_c$ [@Satz:2005hx]. The free energy is a thermodynamical potential of the system which is surrounded by a thermal heat reservoir, such that the temperature is kept fixed by exchanging energy with that reservoir. It measures the amount of available energy of the separated heavy quark anti-quark pair compared to that at different separation distance. On the other hand, the internal energy is a thermodynamical potential of an isolated system so that the energy needed to polarize the heavy quark states should be added. The question, which thermodynamical potential is appropriate as a heavy quark potential at finite temperature, has been controversially discussed in the literature [@Wong:2004zr; @Brambilla:2008cx].

Several years ago, two of the present authors have calculated the strength of the $J/\psi$ wavefunction at the origin as a function of temperature by using a QCD sum rule approach with gluon condensates reliably determined from lattice QCD calculations. The result was then compared with the $J/\psi$ wavefunction obtained by solving the Schrödinger equation either with the free energy or the internal energy as a heavy quark potential [@Lee:2013dca]. It was thus found that the strength of the $J/\psi$ wavefunction as well as the $J/\psi$ mass as a function of temperature obtained from the QCD sum rule calculation follow those obtained from using the free energy as a potential. Meanwhile, the solutions from the internal energy potential significantly overestimate both quantities. This comparison led to the conclusion that the free energy is the appropriate choice for the heavy quark potential at finite temperature rather than the internal energy, which is consistent with recent lattice results extracted from the spectral functions of Wilson line correlators [@Petreczky:2018xuh].

In the potential model, not only the mass of the $J/\psi$, but also that of the $D$ meson is closely related to the heavy quark potential. This can be understood by considering the eigenenergy of an infinitely separated static charm quark pair, which corresponds to twice the $D$ meson mass, as the charmonium will be separated into $D$ and $\bar{D}$ mesons at a large distance with vanishing interaction. The $D$ meson mass at finite temperature hence depends on the behavior of the heavy quark potential at large distance, which is quite different for the free energy and internal energy cases near $T_c$.

In the past, several theoretical studies on the $D$ meson spectrum in a hot or dense medium and its application to heavy-ion collisions have been carried out (see, for example, Refs. [@Hayashigaki:2000es; @Cassing:2000vx; @Mishra:2003se; @Tolos:2007vh; @Suzuki:2015est]). In QCD sum rules, the real part of the current correlator in the $D$ meson channel is related to its imaginary part through a dispersion relation. The former is in the deep Euclidean region computed using the operator product expansion, which leads to an analytic expression with QCD condensates and corresponding Wilson coefficients. The latter imaginary part is expressed as a spectral function composed of physical states with the same quantum number as the D meson. To analyze the sum rules, one then assumes the spectral function to have the simple form of a single ground state and a smooth continuum. The mass of the ground state is then extracted by matching the spectral function to the OPE result. It is thus possible to relate the temperature dependence of the D meson mass to that of the QCD condensates.

The main goal of this paper is to determine what kind of heavy quark potential can model the $D$ meson mass temperature dependence extracted from the QCD sum rule approach. We find that the heavy quark potential at large separation should be composed mostly of the free energy with an addition of a small but non-trivial fraction of the internal energy. Combined with our previous results for the heavy quark potential obtained by comparing the potential model result for the $J/\psi$ to the QCD sum rule calculation, we conclude that the effective heavy quark potential is dominated by the free energy at short distance, while at larger separation, it has a non-negligible fraction of the internal energy.

This paper is organized as follows. In Section \[potentialm\] we introduce the heavy quark potentials from lattice QCD calculations and show the obtained $D$ meson mass as a function of temperature for each of the various potentials. The same $D$ meson mass is calculated from QCD sum rules in Section \[QSR\] and compared to the potential model results in Section \[comparison\]. The conclusion is given in Section \[conclusion\].

$D$ meson mass from heavy quark potential {#potentialm}
=========================================

The strong interaction between the heavy quark and anti-quark is in our approach modeled as a combination of a Coulomb and linear potential. The former is dominant at short, the latter at long distances. The linear potential is furthermore generally responsible for the confinement of hadrons.

![The screening mass $\mu$ scaled by $\sqrt{\sigma}$, extracted from lattice QCD calculations [@Satz:2005hx; @Kaczmarek:1900zz; @Digal:2005ht] (red points) and fitted to Eq.(\[param\]) (blue solid line).[]{data-label="screeningm"}](screeningm.eps){width="8.6"}

The free energy obtained from lattice QCD calculations can be decomposed as,

$$\begin{aligned}
F(r,T)=F_c(r,T)+F_s(r,T),\end{aligned}$$

where $F_c$ and $F_s$, respectively, denote Coulomb-like and the string terms and are expressed as

$$\begin{aligned}
F_c(r,T)&=&-\frac{\alpha}{r}[e^{-\mu r}+\mu r],\nonumber\\
F_s(r,T)&=&\frac{\sigma}{\mu}\bigg[\frac{\Gamma(1/4)}{2^{3/2}\Gamma(3/4)}-\frac{\sqrt{\mu r}}{2^{3/4}\Gamma(3/4)}K_{1/4}[(\mu r)^2]\bigg],\nonumber\\\end{aligned}$$

with $\alpha=\pi/12$, $\sigma=0.445~{\rm GeV^2}$, and $\mu$ being the screening mass [@Satz:2005hx]. In Fig. \[screeningm\] we show $\mu$ extracted from lattice calculations [@Satz:2005hx; @Kaczmarek:1900zz; @Digal:2005ht] as a function of temperature. We parameterize it separately below and above $T_c$ as follows:

$$\begin{aligned}
\frac{\mu}{\sqrt{\sigma}}&=&0.35+0.0034 \exp[(T/T_c)^2/0.22]\nonumber\\
&&~~~~~~~~~{\rm for~} T<T_c,\nonumber\\
\frac{\mu}{\sqrt{\sigma}}&=&0.45+0.5\bigg(\frac{T}{T_c}\bigg) \tanh\big[(T/T_c-0.93)/0.15\big] \nonumber\\
&&~~~~~~~~~{\rm for~} T>T_c.
\label{param}\end{aligned}$$

The two separate parametrizations are smoothly matched at $T_c$, that is, with (almost) the same derivative. Fig. \[screeningm\] shows that the screening mass rapidly increases near $T_c$. If the free energy is interpreted as the heavy quark potential, this sudden increase brings about the lowering of the heavy quark potential, which causes the quarkonia to melt at relatively low temperatures.

The internal energy is derived from thermodynamic relations as $$\begin{aligned}
U(r,T)=F(r,T)+TS(r,T)\nonumber\\
=F(r,T)-T\frac{\partial F(r,T)}{\partial T}.\end{aligned}$$

As mentioned above, the free energy drops sharply near $T_c$ and hence the entropy density defined as $-\partial F/\partial T$ rapidly increases. As a result, if the internal energy is interpreted as the heavy quark potential, a large potential wall is generated near $T_c$, causing the quarkonia to be strongly bound and to dissolve only at higher temperatures.

The heavy quark potential is used in the Schrödinger equation of a charm and anti-charm pair as $$\begin{aligned}
\bigg[2m_c-\frac{1}{m_c}\nabla^2+V(r,T)\bigg]\psi(r,T)=M \psi(r,T),\label{schrodinger1}\end{aligned}$$ where $m_c=1.25$ GeV is the bare charm quark mass and $\psi(r,T)$ the charmonium wave function at temperature $T$. Introducing the potential energy at infinity as $V(r=\infty, T)$, Eq.(\[schrodinger1\]) is modified as [@Karsch:1987pv; @Satz:2005hx] $$\begin{aligned}
\bigg[-\frac{\nabla^2}{m_c}+\widetilde{V}(r,T)\bigg]\psi(r,T)
=-\varepsilon\psi(r,T),
\label{schrodinger2}\end{aligned}$$ where $\widetilde{V}(r,T)\equiv V(r,T)-V(r=\infty, T)$, which vanishes at infinity, while $\varepsilon$ is the $J/\psi$ binding energy at temperature $T$, $$\begin{aligned}
\varepsilon = 2m_c + V(r=\infty,T) - M.\end{aligned}$$

The heavy quark potential in principle also has an imaginary part. Since it however does neither much change the eigenenergy nor the eigenfunction of the Schrödinger equation, we ignore it in the present study [@Laine:2006ns; @Brambilla:2008cx; @Rothkopf:2011db]. At vanishing binding energy, the eigenvalue $M$ of Eq. (\[schrodinger1\]) can be interpreted as the sum of the masses of two open heavy flavors, $D$ and $\bar{D}$ mesons in a hadron gas or dressed charm and anticharm quarks in the QGP,

$$\begin{aligned}
m_D(T)\equiv \frac{M}{2} =m_c+\frac{1}{2}V(r=\infty,T).
\label{dmass}\end{aligned}$$

In other words, if the charm and anticharm quarks that form a quarkonium state are separated from each other by an infinite distance, the energy of the charm pair becomes the energy of two open heavy flavors.

From the above discussion, it is understood that the $D$ meson mass at finite temperature is closely related to the heavy quark potential at infinity. The asymptotic values of the free and internal energies are given by

$$\begin{aligned}
\lim_{r\rightarrow \infty} F(r,T)&=& \frac{\Gamma(1/4)}{2^{3/2}\Gamma(3/4)}\frac{\sigma}{\mu}-\alpha \mu, \nonumber\\
\lim_{r\rightarrow \infty} U(r,T)&=&\frac{\Gamma(1/4)}{2^{3/2}\Gamma(3/4)}\frac{\sigma}{\mu}-\alpha \mu\nonumber\\
&&+ T\frac{d\mu}{dT}\bigg[\frac{\Gamma(1/4)}{2^{3/2}\Gamma(3/4)}\frac{\sigma}{\mu^2}+\alpha\bigg].
\label{eqinfV}\end{aligned}$$

![(a) The free energy and internal energy potentials at infinity and (b) the $D$ meson or dressed charm quark mass for several heavy quark potentials as functions of temperature.[]{data-label="infV"}](infFU.eps){width="8.6"}

![(a) The free energy and internal energy potentials at infinity and (b) the $D$ meson or dressed charm quark mass for several heavy quark potentials as functions of temperature.[]{data-label="infV"}](Dmass.eps){width="8.6"}

The upper panel of Fig. \[infV\] shows the free energy and internal energy potential values at infinity as functions of temperature, obtained from Eqs. (\[param\]) and (\[eqinfV\]). We note that the asymptotic values of the internal energy are very close to those of Refs. [@Satz:2005hx; @Digal:2005ht].

The $D$ meson or dressed charm quark mass is shown in the lower panel of Fig. \[infV\] as a function of temperature for the free energy potential, the internal energy potential, and the averaged potential. One can see in this figure that the behavior of the $D$ meson mass is completely different, depending on the employed potential. If the free energy is used as the potential, the $D$ meson mass decreases with increasing temperature, while the mass sharply increases up to 3.2 GeV near $T_c$ for the internal energy potential . Even if one takes as the heavy quark potential the average of the free and internal energy, the D meson mass increases up to 2.4 GeV around $T_c$.

$D$ meson mass from QCD sum rules {#QSR}
=================================

In this section, we discuss the QCD sum rule approach used to compute the D meson mass at finite temperature. As it is costumary in a QCD sum rule analysis [@Shifman:1978bx; @Shifman:1978by], we start with the correlator of an operator $j(x) = \bar{q}(x) i \gamma_5 c(x)$, which couples to the $D$ meson state of interest, $$\begin{aligned}
\Pi(q_0^2, \bm{q}^2, T) = i \int d^4x e^{iqx} \langle T[ j(x) j^{\dagger}(0) ] \rangle_{T}.
\label{correlator}\end{aligned}$$ In this paper, we will assume that the $D$ is at rest, will hence set $\bm{q} = 0$ and omit this variable for the rest of the discussion for simplicity of notation. The spectral function $\rho(q_0^2, T) = 1/\pi \tanh(q_0/2T) \mathrm{Im} \Pi(q_0^2, T)=1/\pi  \mathrm{Im} \Pi^{\rm R}(q_0^2, T)$, where the subscript R denotes the retarded correlation function, encodes the behavior of the $D$ meson at finite density. With the help of the analyticity of $\Pi^{\rm R}(q_0^2, T)$, the correlator in the deep Euclidean region ($q_0^2 \to - \infty$), calculated through the time ordered correlator, can be related to an integral over $\rho(q_0^2, T)$ via a dispersion relation. Applying furthermore the Borel transform, one obtains the most commonly used version of QCD sum rules, which reads $$\begin{aligned}
\widehat{\Pi}(M^2, T) = \int_{0}^{\infty} ds e^{-s/M^2} \rho(s, T),
\label{sum_rules}\end{aligned}$$ where $\widehat{\Pi}(M^2, T)$ stands for the Borel transformed correlator and $M$ is the so-called Borel mass. For our purposes, we can identify the imaginary part of the time ordered correlator with the spectral density as the thermal factor $\tanh\bigl[\sqrt{s}/(2T)\bigr]$ in the integrand on the right-hand side for the temperature and energy values considered here ($\sqrt{s} \gtrsim m_D$ and $T \lesssim T_c$) is close to 1 and can be safely neglected (see also the discussion in Ref.[@Buchheim:2018kss]).

As mentioned above, the left-hand side of Eq.(\[sum\_rules\]) is computed in the deep Euclidean region. This means that one can rely on the operator product expansion (OPE), which is valid in this regime. We will here not go into the details of this calculation, which have been discussed extensively in the literature (see especially Ref.[@Zschocke:2011aa] for details), but only give the final result to be used in our analysis. The leading order perturbative part of the OPE is given as $$\begin{aligned}
\widehat{\Pi}^{\mathrm{pert}}(M^2, T) = \frac{1}{\pi} \int_{m^2_c}^{\infty} ds e^{-s/M^2} \mathrm{Im} \Pi^{\mathrm{pert}}(s, T)
\label{OPE_pert_1}\end{aligned}$$ with $$\begin{aligned}
\mathrm{Im} \Pi^{\mathrm{pert}}(s, T) = \frac{3}{8 \pi} s \Bigl(1 - \frac{m^2_c}{s}\Bigr)^2  \Bigl(1 + \frac{4}{3} \frac{\alpha_s}{\pi} R(m^2_c/s) \Bigr)
\label{OPE_pert_2}\end{aligned}$$ and $$\begin{aligned}
R(x) &=& \frac{9}{4} + \mathrm{Li}_2(x) + \ln (x) \ln(1-x) \nonumber\\
&&- \frac{3}{2} \ln \Bigl(\frac{1-x}{x}\Bigr) - \ln(1-x) + x \ln \Bigl(\frac{1-x}{x}\Bigr) \nonumber \\
&&- \frac{x}{1-x} \ln (x).
\label{OPE_pert_3}\end{aligned}$$ This term does not depend on the temperature $T$. Next, the non-perturbative scalar condensate terms up to dimension 5 are $$\begin{aligned}
\widehat{\Pi}^{\langle \bar{q} q \rangle}(M^2, T) &=& -m_c \langle \bar{q} q \rangle_{T} e^{-m^2_c/M^2}, \\
\widehat{\Pi}^{\langle \frac{\alpha_s}{\pi} G^2 \rangle}(M^2, T) &=& \frac{1}{12} \Bigl \langle \frac{\alpha_s}{\pi} G^2 \Bigr \rangle_{T} e^{-m^2_c/M^2}, \\
\widehat{\Pi}^{\langle \bar{q} g \sigma G q \rangle}(M^2, T) &=& \frac{1}{2} \Bigl(\frac{m^3_c}{2M^4} - \frac{m_c}{M^2} \Bigr) \nonumber \\
&& \times \langle \bar{q} g \sigma G q \rangle_{T} e^{-m^2_c/M^2},
\label{OPE_non_pert_1}\end{aligned}$$ while the non-perturbative non-scalar condensate terms are $$\begin{aligned}
\widehat{\Pi}^{G_2}(M^2, T) &=& \frac{1}{4}  \Biggl[\Bigl(\frac{7}{6} \ln \frac{\mu^2 m^2_c}{M^4} + 2\gamma_{E} \Bigr) \Bigl(\frac{m^2_c}{M^2} - 1\Bigr) \nonumber \\
&& - 2 \frac{m^2_c}{M^2} \Biggr] G_2(T) e^{-m^2_c/M^2}, \\
\widehat{\Pi}^{F}(M^2, T) &=& \frac{3}{2} \Bigl(\frac{m^2_c}{M^2} - 1\Bigr) F(T) e^{-m^2_c/M^2}, \\
\widehat{\Pi}^{H}(M^2, T) &=&-3 \Bigl( \frac{m^3_c}{2M^4} - \frac{m_c}{M^2} \Bigr) H(T) e^{-m^2_c/M^2}.
\label{OPE_non_pert_2}\end{aligned}$$ The sign in front of the $\gamma_{E}$ (the Euler-Mascheroni constant) term is different from the expression given in Ref.[@Buchheim:2018kss], which can be traced back to a mistake/typo in Eq.(65) of Ref.[@Zschocke:2011aa] and which is corrected here. $G_2(T)$, $F(T)$ and $H(T)$ are related to non-scalar gluon and quark condensates as $$\begin{aligned}
\Bigl \langle \mathcal{ST} \frac{\alpha_s}{\pi}  G^{a \mu}_{\alpha} G^{a \alpha \nu} \Bigr \rangle_T &=& G_2(T) \mathcal{ST}(u_{\mu} u_{\nu}), \nonumber \\
&=& G_2(T)(u_{\mu} u_{\nu} - \frac{1}{4} g_{\mu \nu} ), \\
\langle \mathcal{ST} \bar{q} \gamma_{\mu} i D_{\nu} q \rangle_T  &=& F(T) \mathcal{ST}(u_{\mu} u_{\nu}), \\
\langle \mathcal{ST} \bar{q} i D_{\mu} i D_{\nu} q \rangle_T &=& H(T) \mathcal{ST}(u_{\mu} u_{\nu}),
\label{non_scalar_cond}\end{aligned}$$ where $\mathcal{ST}$ stands for the operation of making the following non-contracted Lorentz indices symmetric and traceless. $u_{\mu}$ is the four-velocity of the heat bath, which we assume to be at rest, hence $u_{\mu} =(1,0,0,0)$. Higher order terms containing dimension 6 four quark condensates were computed in Ref.[@Buchheim:2014rpa], but turned out to be small. We therefore omit them in this work for simplicity. In the QCD sum rule approach considered here, the temperature dependence is assumed to appear only in the QCD condensates and not in their Wilson coefficients. Remembering that the OPE is realized through a division of scales, with low-energy components entering the operators and high-energy contributions expressed as Wilson coefficients (with the boundary roughly at $\Lambda_{QCD}$), it becomes clear that the above assumtion is only valid as long as the temperature is low enough. Usually, one hence regards this OPE to be applicable for temperatures up to and around $T_c$.

For the OPE input parameters at $T = 0$, we employ the values given in Table \[inputparameters\].

                            input parameter                                       value                          reference
  -------------------------------------------------------------------- ---------------------------- ------------------------------------
                                 $m_c$                                     $1.67\,\mathrm{GeV}$             [@Tanabashi:2018oca]
                               $\alpha_s$                                          0.30                     [@Tanabashi:2018oca]
                      $\langle \bar{q}q \rangle_0$                      $(-0.272\,\mathrm{GeV})^3$            [@Aoki:2019cca]
             $\langle \frac{\alpha_s}{\pi} G^2 \rangle_{0}$              $0.012\,\mathrm{GeV}^4$     [@Shifman:1978bx; @Shifman:1978by]
   $\langle \bar{q}g \sigma G q \rangle_0/\langle \bar{q}q \rangle_0$     $0.62\,\mathrm{GeV}^2$             [@Belyaev:1982sa]

  : Vacuum input parameters, given at a renormalization scale of 2 GeV. Note that the $m_c$ used here is the pole mass.[]{data-label="inputparameters"}

To compute the strong coupling constant $\alpha_s$, the expression provided in the PDG [@Tanabashi:2018oca] with $N_f = 4$ and $\Lambda^{(4)}_{\bar{MS}} = 292\,\mathrm{MeV}$ is used. For $\langle \bar{q}g \sigma G q \rangle_0/\langle \bar{q}q \rangle_0$, we use the standard value of $0.8\,\mathrm{GeV}^2$, renormalized at 1 GeV and run it to 2 GeV making use of its anomalous dimension [@Yang:1993bp]. The temperature dependence of the condensates is obtained from different sources. Whenever possible, we have implemented state-of-the-art lattice QCD results. Especially for the most important chiral condensate $\langle \bar{q}q \rangle_T$ (as well as the less important $\langle \frac{\alpha_s}{\pi} G^2 \rangle_{T}$), full QCD, continuum extrapolated results with $2 + 1$ flavors are available at the physical pion mass. Specifically, we make use of the lattice data given in Ref.[@Borsanyi:2010bp] for $\langle \bar{q}q \rangle_T$ and those of Ref.[@Bazavov:2014pvz] for $\langle \frac{\alpha_s}{\pi} G^2 \rangle_{T}$, according to the prescription discussed in Ref.[@Gubler:2018ctz]. For $\langle \bar{q}g \sigma G q \rangle_T$, which so far has not been studied extensively on the lattice, we assume $\langle \bar{q}g \sigma G q \rangle_T/\langle \bar{q}q \rangle_T$ to be constant and can therefore extract the $\langle \bar{q}g \sigma G q \rangle_T/\langle \bar{q}q \rangle_T$ temperature dependence directly from $\langle \bar{q}q \rangle_T$. The above assumption is supported by the lattice QCD calculation of Ref.[@Doi:2004jp]. The non-scalar quark condensates encoded in $G_2(T)$, $F(T)$ and $H(T)$ are the least well known input parameters for the sum rules. To estimate them, we apply the hadron resonance gas model, which assumes that the effect of finite temperature can be described by a gas of non-interacting pions (and other light hadrons). This assumtion should be approximately valid below $T_c$, where hadronic degrees of freedom dominate, but becomes questionable around or above $T_c$, where quarks and gluons start to appear. For concrete formulas and other details, we refer the interested reader to Ref.[@Gubler:2018ctz], and just mention here that we have included all small-mass pesudo-scalar particles, namely pions, kaons and the eta in the actual claculation. To obtain $G_2(T)$, one furthermore needs the value of the strong coupling constant at finite temperature, $\alpha_s(T)$. To extract it, we use the two-loop perturbative running coupling, as given in Ref.[@Morita:2007hv] (see also Ref.[@Kaczmarek:2004gv]).

With all the OPE input determined, we next discuss how to analyze the sum rules of Eq.(\[sum\_rules\]). As it is standard, we assume the spectral function to have the so-called “pole + continuum” form (see, however, Ref.[@Gubler:2010cf] for an alternative approach), $$\begin{aligned}
\rho(s, T) &=&
\lambda_D \delta(s - m^2_D) \nonumber \\
&& + \theta(s - s_{th}) \frac{1}{\pi} \mathrm{Im} \Pi^{\mathrm{pert}}(s, T),
\label{pole_cont}\end{aligned}$$ which can be expected to be justified when the $D$ meson ground state peak dominates the spectral function at low energy and the perturbative expression $\Pi^{\mathrm{pert}}(s, T)$ descibes the spectral function above the threshold parameter $s_{th}$ reasonably well. Some care is however needed when assessing the validity of the above assumption. While it may be sufficiently accurate in vacuum ($T = 0$), it will eventually break down once the peak starts to dissolve into the continuum at some temperature above $T_c$. Analyses based on Eq.(\[pole\_cont\]) hence can only be trusted for temperatures up to and around $T_c$.

Sustituting Eq.(\[pole\_cont\]) into Eq.(\[sum\_rules\]), we have $$\begin{aligned}
\widehat{\Pi}(M^2, T) &=& \lambda_D e^{-m^2_D/M^2} \nonumber \\
&& + \frac{1}{\pi} \int_{s_{th}}^{\infty} ds e^{-s/M^2} \mathrm{Im} \Pi^{\mathrm{pert}}(s, T).
\label{sum_rules_2}\end{aligned}$$ Therefore, we define $$\begin{aligned}
&& \widetilde{\Pi}(M^2, T, s_{th}) \equiv \nonumber \\
&& \widehat{\Pi}(M^2, T) - \frac{1}{\pi} \int_{s_{th}}^{\infty} ds e^{-s/M^2} \mathrm{Im} \Pi^{\mathrm{pert}}(s, T),
\label{sum_rules_3}\end{aligned}$$ and can compute the $D$ meson mass as $$\begin{aligned}
m_D(M^2, T, s_{th}) = \sqrt{-\frac{ \partial \widetilde{\Pi}(M^2, T, s_{th})/\partial(1/M^2)}{\widetilde{\Pi}(M^2, T, s_{th})}},
\label{sum_rules_4}\end{aligned}$$ and subsequently the residue $\lambda_D$ as $$\begin{aligned}
\lambda_D (M^2, T, s_{th}) = \widetilde{\Pi}(M^2, T, s_{th}) e^{m^2_D(M^2, T, s_{th})/M^2}.
\label{sum_rules_5}\end{aligned}$$ We, however, still need to fix $M$ and $s_{th}$, which are artificial parameters of the sum rule approach. Starting with $M$, we define a so-called Borel window, for which the approximations used in the sum rule analysis can be expected to be valid. The lower boundary of the Borel window, $M_{\mathrm{min}}$, is determined from the condition of sufficient convergence of the OPE, namely $$\begin{aligned}
\frac{\widehat{\Pi}^{\mathrm{OPE}}_{\mathrm{dim\,5\,terms}}(M^2, T)}{ \widehat{\Pi}^{\mathrm{OPE}}_{\mathrm{all\,terms}}(M^2, T) } < 0.1.
\label{Borel_window_1}\end{aligned}$$ Here, $\widehat{\Pi}^{\mathrm{OPE}}_{\mathrm{dim\,5\,terms}}(M^2, T)$ is the sum of the terms of condensates with the highest dimension considered in this work, these are $\widehat{\Pi}^{\langle \bar{q} g \sigma G q \rangle}(M^2, T)$ and $\widehat{\Pi}^{e_2}(M^2, T)$. For the upper boundary, $M_{\mathrm{max}}$, we employ the condition that the pole contribution in the integral of Eq.(\[sum\_rules\]) should dominate the sum rule. Specifically, we have $$\begin{aligned}
\frac{  \int_{0}^{s_{th}} ds e^{-s/M^2} \rho(s, T)  }{ \int_{0}^{\infty} ds e^{-s/M^2} \rho(s, T) }
= \frac{ \widetilde{\Pi}(M^2, T, s_{th}) }{ \widehat{\Pi}(M^2, T) } > 0.5.
\label{Borel_window_2}\end{aligned}$$ With the range of $M$ determined as $M_{\mathrm{min}} < M < M_{\mathrm{max}}$, one can compute averages of $m_D(M^2, T, s_{th})$ and $\lambda_D (M^2, T, s_{th})$. We have $$\begin{aligned}
\overline{m}_D(T, s_{th}) &=& \frac{1}{M_{\mathrm{max}} - M_{\mathrm{min}}} \nonumber \\
&& \times \int_{M_{\mathrm{min}}}^{M_{\mathrm{max}}} dM m_D(M^2, T, s_{th}), \\
\overline{\lambda}_D (T, s_{th}) &=& \frac{1}{M_{\mathrm{max}} - M_{\mathrm{min}}} \nonumber \\
&& \times \int_{M_{\mathrm{min}}}^{M_{\mathrm{max}}} dM \lambda_D (M^2, T, s_{th}).
\label{averages}\end{aligned}$$ Finally, $s_{th}$ is determined as $s_{th}^0$, which minimizes the function $$\begin{aligned}
&& d(T, s_{th}) =  \frac{1}{M_{\mathrm{max}} - M_{\mathrm{min}}} \nonumber \\
&& \times \int_{M_{\mathrm{min}}}^{M_{\mathrm{max}}} dM \Bigl[ m_D(M^2, T, s_{th}) - \overline{m}_D(T, s_{th}) \Bigr]^2.
\label{determine_s_th}\end{aligned}$$ This ensures that, within the Borel window, the dependence of $m_D(M^2, T, s_{th})$ on $M$ is as small as possible. For illustration, we show in Fig.\[Borel\_curve\_comp\] $m_D(M^2, T, s_{th}^0)$ for $T = 0$ and $T = 150$ MeV as a funtion of $M$. The respective positions of $M_{\mathrm{min}}$ and $M_{\mathrm{max}}$ are indicated as vertical arrows.

![The $D$ meson masses computed from Eq.(\[sum\_rules\_4\]) at $s_{th}^0$, which gives the maximum flatness, in the vacuum and at $T = 150$ MeV. The vertical arrows indicate the locations of the Borel window boundaries, $M_{\mathrm{min}}$ and $M_{\mathrm{max}}$.[]{data-label="Borel_curve_comp"}](mD_Mdep_comp_2.eps){width="8.6"}

We have carried out two different analyses, one in which we let the Borel window depend on $T$ (shown in Fig.\[Borel\_curve\_comp\]) and another one in which we keep it fixed at its vacuum location.

![The values of $\overline{m}_D(T, s_{th}^0)$ and $\overline{\lambda}_D (T, s_{th}^0)$ as a function of temperature $T$. The red curve (horizontal crosses) shows the analysis results when the Borel window is determined at $T=0$ and kept fixed with increasing temperature. The blue curve (diagonal crosses) gives the results for a temperature dependent Borel window.[]{data-label="Borel_window_comp"}](mD_T_comp_2.eps){width="8.6"}

![The values of $\overline{m}_D(T, s_{th}^0)$ and $\overline{\lambda}_D (T, s_{th}^0)$ as a function of temperature $T$. The red curve (horizontal crosses) shows the analysis results when the Borel window is determined at $T=0$ and kept fixed with increasing temperature. The blue curve (diagonal crosses) gives the results for a temperature dependent Borel window.[]{data-label="Borel_window_comp"}](res_T_comp_2.eps){width="8.6"}

The two versions are shown in Fig.\[Borel\_window\_comp\], where it is seen that the $D$ meson mass increases or remains roughly constant below and around $T_c \simeq 155$ MeV [@Borsanyi:2010bp] and starts to drop for higher temperatures. The residue $\lambda_{D}$, on the other hand, decreases monotonically with increasing temperature and approaches zero towards $T = 180$ MeV. Comparing the two analyses, it is observed that only the $T$ dependent Borel window case leads to an increase of the $D$ meson mass. It also gives a slightly slower decrease of the residue $\lambda_D$.

It is worth thinking about a physically intuitive picture for the above results. First, let us consider the potential increase of the $D$ meson mass below and around $T_c$, which seems to be consistent with the reconstructed spectral functions given in the recent lattice QCD study of Ref.[@Kelly:2018hsi], the lattice data quality however not being good enough to draw any definite conclusions. The most natural explanation for this behavior can be given by the quark model picture proposed in Ref.[@Park:2016xrw], where the increase of the $D$ meson mass was related to the decreasing constituent quark mass caused by the partial restoration of chiral symmetry. Due to the decreasing quark mass, its wave function spreads out to larger distances and therefore feels the linear confining potential, leading to an averall increase of the $D$ meson mass. As in the present study the chiral symmetry gets partially restored due to the decrease of the chiral condensate, the same quark model picture provides a natural explanation for an increasing $D$ meson mass behavior around $T_c$. Next, we examine the decreasing $D$ meson mass at temperatures around and above $T_c$. To consider this behavior, one should take into account the decreasing residue value $\lambda_{D}$, shown in the lower plot of Fig.\[Borel\_window\_comp\]. This decrease means that the contribution of the ground state peak to the spectral function is reduced, which can be interpreted as the gradual melting of the $D$ meson into the continuum, which is expected to happen at some (unknown) temperature above $T_c$. Therefore, if the spectral function at $T > T_c$ is gradually changed from a narrow peak to a continuum, the “pole + continuum” assumption of Eq.(\[pole\_cont\]) breaks down and the notion of a $D$ meson mass looses any meaning. It should however be kept in mind that both the division of scales of the OPE and the estimation of some of the non-scalar QCD condensates become questionable above $T_c$. The physical interpretation of the drop shown in the upper plot of Fig.\[Borel\_window\_comp\] for $T > T_c$ hence is not completely clear because of the limitations of the QCD sum rule method at such temperatures, but is likely related to the melting of the $D$ meson state into the continuum. In the next section, we will discuss further physical interpretations of the QCD sum rule results in terms of the quark model approach considered in the previous section.

Matching the heavy quark potential to QCD sum rules {#comparison}
===================================================

![The $D$ meson or dressed charm quark mass for several heavy quark potentials as functions of temperature in comparison with the QCD sum rule results. $T_c$ is taken to be 155 MeV. In the potential model $m_c$ is set to 1.29 GeV. For the results from Ref. [@Montana:2020lfi], effects from spectral widths are shown as error bars. All results are shifted to have the common physical $D$ meson mass in vacuum.[]{data-label="philipp-juan"}](philipp.eps){width="8.6"}

Let us now study the consequences of requiring that the temperature dependence of the D meson mass obtained from the potential model given in Eq. (\[dmass\]) is consistent with that from QCD sum rules. As can be seen in Fig.\[philipp-juan\], assuming a constant or temperature dependent Borel window requires that the heavy quark potential should have a nontrivial contribution from the internal energy in addition to the dominant free energy part. The discrepancy above $T_c$ should not be taken too seriously as the $D$ meson mass calculations are not reliable at these temperatures, as mentioned in the previous section. Also shown in the figure is the recent results from Ref.[@Montana:2020lfi] where the $D$ meson mass was calculated at finite temperature in the imaginary-time formalism by using a hadronic effective field theory with chiral and heavy-quark spin symmetries. Although the central value for the $D$ meson mass slightly decreases as the temperature approaches the critical temperature, the difference compared with the QCD sum rule result does not qualitatively change our discussion of the heavy quark potential of this section. We find that fitting the results in the potential model to the $D$ meson mass obtained from QCD sum rules with either a constant or temperature dependent Borel window requires the heavy quark potential to be composed of 18% and 21% of internal energy potentials, respectively, with the remaining contribution coming from the free energy potential. The results from the hadronic effective model are fitted with 15% of the internal energy, if one fits the potential to the central pole mass.

On the other hand, a previous determination of the heavy quark potential obtained by comparing the potential model result for the $J/\psi$ to a QCD sum rule calculation, led to the conclusion that the heavy quark potential should be dominated by the free energy [@Lee:2013dca]. The two findings are, however, not necessarily inconsistent, as combining them, one can naturally conclude that the effective heavy quark potential is dominated by the free energy at short distance, while at larger separation, it will have a non-trivial fraction of internal energy. This is in line with the interpretation given in Ref. [@Satz:2015jsa]. For short distances, the interquark potential should be affected by the gluon dynamics exchanged between the two heavy quarks, which is encoded in the free energy. On the other hand, as the separation between the heavy quarks becomes large, the polarization of the heavy quark itself becomes important, which is related to the D meson, but at the same time introduces additional entropy contribution included in the internal energy. The same effect can be simulated with using just the free energy potential but a larger effective heavy quark mass at larger distance. The introduction of additional internal energy contribution to the potential energy will become important when one calculates the binding energy and wave function of the excited states which are of larger size compared to the ground state.

![(a) Three different heavy quark potential energies, (b) the residue and (c) the mass of $J/\psi$ as a function of temperature for each heavy quark potential energy, compared with the results from QCD sum rules [@Lee:2013dca].[]{data-label="jpsi"}](transitional-potential.eps){width="8.6"}

![(a) Three different heavy quark potential energies, (b) the residue and (c) the mass of $J/\psi$ as a function of temperature for each heavy quark potential energy, compared with the results from QCD sum rules [@Lee:2013dca].[]{data-label="jpsi"}](residue-compare.eps){width="8.6"}

![(a) Three different heavy quark potential energies, (b) the residue and (c) the mass of $J/\psi$ as a function of temperature for each heavy quark potential energy, compared with the results from QCD sum rules [@Lee:2013dca].[]{data-label="jpsi"}](mass-compare.eps){width="8.6"}

As an illustration we introduce a heavy quark potential, which starts with the free energy at short distance and then smoothly transits to a combinational potential energy composed of 80% of the free energy and 20% of the internal energy with increasing interquark distance, as shown in Fig. \[jpsi\](a). In Fig. \[jpsi\](b) and (c), we show the result for the wavefunction at the origin and the $J/\psi$ mass obtained by solving the Schrödinger equation of Eq.(\[schrodinger2\]) by using three different potentials: the free energy potential, the combinational potential, and the transitional potential. Also shown in the figures are the results obtained from QCD sum rules [@Lee:2013dca]. The results for the free energy are slightly different from those reported in Ref. [@Lee:2013dca] because of the parametrization of Eq. (\[param\]). First, one notes that results obtained by using the combinational potential energy with 20% of the internal energy deviates from the temperature dependence obtained from the QCD sum rule approach. However, comparing the results between the free energy and the transitional potential, shows that the residue and mass of the $J/\psi$ from the latter better reproduces the QCD sum rule results at least up to the critical temperature where the $D$ meson mass is investigated in this study. Furthermore, the quoted QCD sum rule result is obtained within the zero width approximation. If one allows for the increase of thermal width in the assumed spectral density the obtained mass shift can be expected to be smaller to satisfy the sum rule [@Morita:2007hv]. Hence, the sum rule calculation should be considered as an upper bound for the mass shift as well as the modification of the residue. As the width will rapidly increase above $T_c$, the agreement between the QCD sum rule result and the potential model result with the mixed potential is expected to persist even at these temperatures. Our comparison therefore suggests the heavy quark potential to be composed purely of the free energy at short distance while at larger separation it has a non-trivial contribution from the internal energy. This configuration can reproduce both the $J/\psi$ and the $D$ meson masses at finite temperature obtained from QCD sum rules up to and around the critical temperature.

Summary and Conclusions {#conclusion}
=======================

We have in this paper studied the heavy quark - anti-quark potential at finite temperature, which is crucial for understanding the quarkonium production in relativistic heavy-ion collisions, through the behavior of the $D$ meson mass obtained by using QCD sum rules. Based on the observation that the total energy of a heavy quark bound state pair at infinity can be interpreted as twice the D meson mass at finite temperature, it becomes possible to compute the heavy quark potential at a large distance from the thermal behavior of D mesons. Our QCD sum rule results for the D meson at finite temperature suggest that the heavy quark potential at large distance should be composed dominantly of the free energy but with a contribution of the internal energy of about 20%. Combined with a similar comparison for the $J/\psi$ channel at finite temperature, we conclude that the heavy quark potential should have a different fraction of internal energy as a function of the interquark distance. Summarizing the result, we find that to reproduce the behavior of both the $J/\psi$ and $D$ meson from QCD sum rules without contradiction at least up to and around $T_c$, the heavy quark potential should be composed of the free energy at short distance, while attaining a non-trivial fraction of the internal energy at larger distances so that at infinity the fraction becomes about 20%.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors acknowledge useful discussions with E. Bratkovskaya and J. M. Torres-Rincon. This work was supported by by Samsung Science and Technology Foundation under Project Number SSTF-BA1901-04. P.G. is supported by Grant-in-Aid for Early-Carrier Scientists No. JP18K13542 and the Leading Initiative for Excellent Young Researchers (LEADER) of the Japan Society for the Promotion of Science (JSPS).

[99]{}

T. Matsui and H. Satz, Phys. Lett.  B [**178**]{}, 416 (1986). A. Rothkopf, arXiv:1912.02253 \[hep-ph\]. H. Satz, Eur. Phys. J. C [**75**]{}, no. 5, 193 (2015).

S. H. Lee, K. Morita, T. Song and C. M. Ko, Phys. Rev. D [**89**]{}, no. 9, 094015 (2014). O. Kaczmarek, F. Karsch, P. Petreczky and F. Zantow, Phys. Lett. B [**543**]{}, 41 (2002). H. Satz, J. Phys. G [**32**]{}, R25 (2006). C. Y. Wong, Phys. Rev. C [**72**]{}, 034906 (2005).

N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, Phys. Rev. D [**78**]{}, 014017 (2008).

P. Petreczky, A. Rothkopf and J. Weber, Nucl. Phys. A **982**, 735-738 (2019).

A. Hayashigaki, Phys. Lett. B **487**, 96-103 (2000). W. Cassing, E. Bratkovskaya and A. Sibirtsev, Nucl. Phys. A **691**, 753-778 (2001). A. Mishra, E. Bratkovskaya, J. Schaffner-Bielich, S. Schramm and H. Stoecker, Phys. Rev. C **69**, 015202 (2004). L. Tolos, A. Ramos and T. Mizutani, Phys. Rev. C [**77**]{}, 015207 (2008). K. Suzuki, P. Gubler and M. Oka, Phys. Rev. C [**93**]{}, no. 4, 045209 (2016). O. Kaczmarek, Eur. Phys. J. C [**61**]{}, 811 (2009). S. Digal, O. Kaczmarek, F. Karsch and H. Satz, Eur. Phys. J. C [**43**]{}, 71 (2005). F. Karsch [*et al.*]{}, Z. Phys.  C [**37**]{}, 617 (1988). M. Laine, O. Philipsen, P. Romatschke and M. Tassler, JHEP [**0703**]{}, 054 (2007). A. Rothkopf, T. Hatsuda and S. Sasaki, Phys. Rev. Lett.  [**108**]{}, 162001 (2012).

M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B [**147**]{}, 385 (1979). M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B [**147**]{}, 448 (1979). T. Buchheim, T. Hilger, B. Kämpfer and S. Leupold, J. Phys. G [**45**]{}, no. 8, 085104 (2018). S. Zschocke, T. Hilger and B. Kämpfer, Eur. Phys. J. A [**47**]{}, 151 (2011). T. Buchheim, T. Hilger and B. Kämpfer, Phys. Rev. C [**91**]{}, 015205 (2015). M. Tanabashi [*et al.*]{} \[Particle Data Group\], Phys. Rev. D [**98**]{}, no. 3, 030001 (2018). K. C. Yang, W. Y. P. Hwang, E. M. Henley and L. S. Kisslinger, Phys. Rev. D [**47**]{}, 3001 (1993). S. Aoki [*et al.*]{} \[Flavour Lattice Averaging Group\], arXiv:1902.08191 \[hep-lat\]. V. M. Belyaev and B. L. Ioffe, Sov. Phys. JETP [**56**]{}, 493 (1982) \[Zh. Eksp. Teor. Fiz.  [**83**]{}, 876 (1982)\]. S. Borsanyi [*et al.*]{} \[Wuppertal-Budapest Collaboration\], JHEP [**1009**]{}, 073 (2010). A. Bazavov [*et al.*]{} \[HotQCD Collaboration\], Phys. Rev. D [**90**]{}, 094503 (2014). P. Gubler and D. Satow, Prog. Part. Nucl. Phys.  [**106**]{}, 1 (2019). T. Doi, N. Ishii, M. Oka and H. Suganuma, Phys. Rev. D [**70**]{}, 034510 (2004). K. Morita and S. H. Lee, Phys. Rev. C [**77**]{}, 064904 (2008). O. Kaczmarek, F. Karsch, F. Zantow and P. Petreczky, Phys. Rev. D [**70**]{}, 074505 (2004) Erratum: \[Phys. Rev. D [**72**]{}, 059903 (2005)\]. P. Gubler and M. Oka, Prog. Theor. Phys.  [**124**]{}, 995 (2010). A. Kelly, A. Rothkopf and J. I. Skullerud, Phys. Rev. D [**97**]{}, no. 11, 114509 (2018). A. Park, P. Gubler, M. Harada, S. H. Lee, C. Nonaka and W. Park, Phys. Rev. D [**93**]{}, no. 5, 054035 (2016). G. Montaña, À. Ramos, L. Tolós and J. M. Torres-Rincon, arXiv:2001.11877 \[hep-ph\].
