1 Introduction

Heavy quark-antiquark bound states, quarkonia, like charmonium and bottomonium, serve as a good thermometer for the quark-gluon-plasma (QGP) created in ultra-relativistic heavy ion (AA) collisions [1]. This is due to the fact that they are produced during the early stages of the collisions and participate the entire evolution of the QGP. Some of them will dissociate into heavy quarks when experiencing temperature increasing while some can remain as bound states due to their hierarchically small sizes and large binding energies [2, 3]. How to theoretically understand the suppression of quarkonia yields observed in AA collisions compared to that in the proton-proton (pp) collisions [4,5,6], however, remains unclear. This is mainly due to the complicated entanglement between the cold and hot nuclear effects [7]. Given that the suppression resides in the non-perturbative regime, naturally we want to interpret it using lattice QCD calculations. Besides, the experiments at RHIC and LHC reveal that open heavy mesons show an unexpectedly substantial elliptic flow, comparable to that of light mesons [8,9,10]. This indicates that heavy quarks flow as efficiently as light quarks do [11,12,13]. Trying to interpret this requires a modeling of the heavy quark diffusion in a hot and dense medium, which is challenging in perturbation theory [14,15,16,17], but feasible from lattice QCD. Studies carried out in the heavy quark mass limit can be found in [18,19,20,21,22,23] in both quenched approximation and (2+1)-flavor lattice QCD. Accessing the first mass suppressed correction to the diffusion coefficient are also attempted [22, 24,25,26].

Tackling the above puzzles relies on the knowledge about the spectral function of quarkonia, which can be extracted from the fully relativistic Euclidean mesonic correlation functions. The spectral function in the vector channel contains all information about the in-medium hadron properties like the dissociation temperatures and heavy quark diffusion, and was intensely studied in the past decade, see [27,28,29,30,31,32,33] for a selection. While the one in the pseudoscalar channel provides an easier access to the fate of the bound states due to the vanishing of the transport peak [34, 35]. This was investigated in the quenched case in [36]. A review of current status of lattice studies on heavy quarks and quarkonium in extreme conditions can be found in [37,38,39].

In this work we aim to provide a (2+1)-flavor lattice QCD calculation for the charmonium and bottomonium correlation functions in pseudoscalar channel, \(G_{\textrm{PS}}(\tau )\), at two temperatures T=220 MeV and T=251 MeV. The Euclidean correlator is related to the spectral function via following equation

$$\begin{aligned} G_{\textrm{PS}}(\tau ) \; = \; \int _0^\infty \frac{\textrm{d}\omega }{\pi } \rho _{\textrm{PS}} (\omega ) \frac{\cosh \left( \left( \frac{1}{2T} - \tau \right) \omega \right) }{\sinh \left( \frac{\omega }{2 T} \right) }\;. \end{aligned}$$
(1)

Due to statistical errors in \(G_{\textrm{PS}}(\tau )\) and a finite resolution of the lattice, the solution of \(\rho _{\textrm{PS}}\) is not unique. However, the spectral reconstruction can be constrained with further inputs. In this study we model the spectral function with a perturbatively inspired ansatz, interpolating between the thermal part around the threshold obtained from pNRQCD [40] and the vacuum part well above the threshold [36].

This article is organized as follows. We begin with the lattice setup and clover mass tuning, which are outlined in Sect. 2. Section 3 presents the numerical data of lattice correlators, while Sect. 4 details the process of obtaining perturbative spectral functions. In addition, Sect. 5 provides an explanation of the spectral reconstruction procedure. Finally, we summarize our findings and present future outlook in the concluding section.

2 Lattice Parameters and Mass Tuning

Quarkonium correlators are obtained at various temperatures, as specified in Table 1, using clover-improved Wilson fermions as valence quarks. The lattice correlation functions are computed on large and fine gauge field configurations generated with \(N_f=2+1\) and with the Highly Improved Staggered Quark (HISQ) action [41] and tree-level improved Lüscher-Weisz gauge action [42, 43]. The sea quark masses are tuned to the physical strange quark mass, \(m_s\), and degenerate up, down quark masses \(m_l=m_s/5\), corresponding to \(m_\pi \simeq 320\) MeV. The fermion valence action is tadpole-improved by incorporating clover-improvement, which reduces the \({\mathcal {O}}(a)\) cut-off effects. To this aim, the Sheikholeslami-Wohlert coefficient [44] is chosen as \(c_{SW}=1/u^3_0\), where \(u_0\) is the tadpole factor determined by taking the fourth root of the plaquette expectation value calculated on our HISQ lattices.

Table 1 Lattice parameters for \(N_f\)=2+1, \(m_l=m_s/5\) HISQ configurations for different temperatures
Full size table

The masses of the lattice calculated correlators are determined by the hopping parameter \(\kappa \), which should be tuned to match the quarkonium meson mass spectrum observed in experiments [46]. In the heavy quark sector, the tuning is performed with respect to the spin-averaged charmonium \(m_{c\bar{c}}=(m_{\eta _{c}}+3m_{J/\psi })/4\) and bottomonium \(m_{b\bar{b}}=(m_{\eta _{b}}+3m_{\Upsilon })/4\). Figure 1 illustrates the charm quark mass tuning, where \(am^{phy}_{c\bar{c}}\) represents the physical mass in lattice units. The values of \(\kappa \) that are tuned for charm and bottom are 0.13164 and 0.11684, respectively.

Fig. 1
figure 1

Charm quark mass tuning on the mixed action. The solid line represents the experimental value of the spin-averaged charmonium in lattice units taken from the PDG [46]. Whereas triangles represent masses from lattice calculations for different values of \(\kappa \) in lattice units

Full size image

3 Correlators

We present the pseudoscalar channel charmonium and bottomonium correlators in Fig. 2. Instead of the correlator itself, we used the ratio \(G_{\textrm{PS}}(\tau )/G^{\textrm{free}}_{\textrm{PS}}(\tau )\) for a better visibility of the behavior as a function of \(\tau \) at various temperatures, where \(G^{\textrm{free}}_{\textrm{PS}}(\tau )\) reads [36]

$$\begin{aligned} \frac{ G_{\textrm{PS}}^{\textrm{free}}(\tau ) }{ m^2(\bar{\mu }^{\textrm{ref}}) } \equiv \int _{2 M_q}^{\infty } \frac{\textrm{d}\omega }{\pi } \, \biggl \{\frac{N_{\textrm{c}}\omega ^2}{8\pi } \tanh \Bigl ( \frac{\omega }{4T} \Bigr ) \sqrt{1 - \frac{4 M_q^2}{\omega ^2}}\biggr \}\, \frac{\cosh \left( \left( \frac{1}{2T} - \tau \right) \omega \right) }{\sinh \left( \frac{\omega }{2T}\right) }, \end{aligned}$$
(2)

where \(m(\bar{\mu }^{\textrm{ref}})\) denotes the running mass at reference scale \(\bar{\mu }^{\textrm{ref}}\)=2 GeV. We choose the values of \(M_q\) to be 1.5 GeV for charmonium and 4.7 GeV for bottomonium. One can see that the charmonium suffers more thermal effects than the bottomonium does, because the charmonium correlators at high temperatures deviate more from the low-temperature ones.

Fig. 2
figure 2

The ratio of pseudoscalar lattice correlator with its corresponding free correlator at different temperatures for charmonium (left) and bottomonium (right). Statistical errors are estimated using the Jackknife procedure

Full size image

4 Perturbative Spectral Functions in PS-Channel

To construct the quarkonium spectral functions in the pseudoscalar channel we follow the strategy developed in the quenched approximation [36]. For \(\omega \gg 2 M\), the spectral function is insensitive to temperature effects, therefore it is replaced by its vacuum counterpart. The vacuum spectral function, normalised by \(\omega ^2 m^2(\bar{\mu })\), in the pseudoscalar channel reads

$$\begin{aligned} \left. \frac{ \rho _{\textrm{PS}}(\omega ) }{\omega ^2 m^2(\bar{\mu })} \right| ^{\textrm{vac}} \; \equiv \; \frac{N_{\textrm{c}}}{8\pi }\, \tilde{R}_{\textrm{c}}^{p}(\omega , \bar{\mu })\;. \end{aligned}$$
(3)

Here \(\tilde{R}\)-function is valid up to \({\mathcal {O}}(\alpha _{\textrm{s}}^3)\) and it is defined in Ref. [36]. One can also express the spectral function \(\rho ^{}_{\textrm{PS}}\) normalized by the pole mass M as follows:

$$\begin{aligned} \left. \frac{ \rho _{\textrm{PS}}(\omega ) }{\omega ^2 M^2} \right| ^{\textrm{vac}} \; \equiv \; \frac{N_{\textrm{c}}}{8\pi }\, {R}_{\textrm{c}}^{p}(\omega ) \;. \end{aligned}$$
(4)

where R is again given in Ref. [36]. However, this expansion does not converge as the correction term becomes increasingly large as \(\omega \) increases. In contrast, the form in Eq. (3), obtained by re-expanding the pole mass in terms of the \(\overline{MS}\) mass in Eq. (4), leads to smaller corrections at large \(\omega \). Following Ref. [36] we can then estimate the pole mass self-consistently by equating Eq. (3) and Eq. (4), which is given by,

$$\begin{aligned} M_x^2 \equiv m^2(\bar{\mu })\, \left. \frac{\tilde{R}^p_{\textrm{c}}(\omega ,\bar{\mu })}{R^p_{\textrm{c}}(\omega )} \right| _{\bar{\mu }=\omega ,~ \omega = x M_x}\;. \end{aligned}$$
(5)

Here the running of the \(\overline{MS}\) mass has been calculated using 5-loop running, starting from \(m(\bar{\mu }= 2 \, \text {GeV})\). Table 2 shows the values of \(M_x\) obtained at \({\mathcal {O}}(\alpha _{\textrm{s}}^2)\).

Table 2 The quark masses relevant for the perturbative spectral functions for full QCD, utilizing 5-loop running for \(\alpha _{\textrm{s}}\) and setting \({\Lambda _{{\overline{\textrm{MS}}}}}\) to 0.339 GeV [47]
Full size table

Now we consider thermal part of the spectral function around the threshold i.e. \(\omega \approx 2 M\). In the heavy quark limit, relativistic effects in quark-antiquark bound states become negligible, allowing for a non-relativistic potential description. For separations much smaller than 1/gT, the potential is effectively a zero temperature potential. However, as the separation approaches \(\sim 1/gT\), the contribution of soft gluon exchange becomes significant, requiring a Hard Thermal Loop resummation. As a result, the potential in this region becomes temperature dependent and acquires a complex component. In this non-relativistic regime, one can get the pseudoscalar spectral function from the vector channel as

$$\begin{aligned} \rho ^{\textrm{pNRQCD}}_{\textrm{PS}} =\frac{M^2}{3} \rho ^{\textrm{pNRQCD}}_{\textrm{V}}\;. \end{aligned}$$
(6)

The vector channel spectral function is given by

$$\begin{aligned} \rho ^{\textrm{pNRQCD}}_{\textrm{V}}(\omega ) =\frac{1}{2} \Bigl ( 1 - e^{-\frac{\omega }{T }}\Bigr ) \int _{-\infty }^{\infty } \textrm{d} t \, e^{i \omega t} \; C_{>}(t;{\vec {0},\vec {0}})\;, \end{aligned}$$
(7)

where \(C_{>}\) is a Wightman function and it is obtained by solving the following Schrödinger equation

$$\begin{aligned} \biggl \{ i \partial _t - \biggl [ 2 M + V_{\textrm{T}}(r) - \frac{\nabla ^2_{\vec {r}}}{M} \biggr ] \biggr \} \, C_{>}^V(t;{\vec {r},\vec {r'}}) = 0 \;, \quad t\ne 0 \;, \end{aligned}$$
(8)

with initial conditions \(C_{>}^V(0;{\vec {r},\vec {r'}}) = 6 N_{\textrm{c}}\, \delta ^{(3)}({\vec {r}-\vec {r'}})\;\). The potential \(V_{\textrm{T}}(r)\) for positive t has the following form [48,49,50]

$$\begin{aligned} V_{\textrm{T}}(r) = -\alpha _{\textrm{s}}C_F \biggl [ m_{\textrm{D}}+ \frac{\exp (-m_{\textrm{D}}r)}{r} \biggr ] - {i \alpha _{\textrm{s}}C_{\textrm{F}}T } \, \phi (m_{\textrm{D}}r) + {\mathcal {O}}(\alpha _{\textrm{s}}^2)\;. \end{aligned}$$
(9)

Here \(\phi (x)\) is defined to be

$$\begin{aligned} \phi (x) \equiv 2 \int _0^\infty \frac{\textrm{d} z \, z}{(z^2 +1)^2} \biggl [1 - \frac{\sin (z x)}{zx}\biggr ]. \end{aligned}$$
(10)

Here \(\alpha _{\textrm{S}}\) is the thermal coupling and \(m_{{\textrm{D}}}\) \((\sim gT)\) is the Debye mass. Physically, the real part of the potential is related to Debye screening in the plasma. It affects the thermal mass shift of bound states. On the other hand, the imaginary part originates from Landau damping of space like gauge fields, which causes the bound state peak to broaden. At a much lower \(\omega \) than the threshold, the Schrödinger description described above overestimates the spectral function, as this formalism is no longer valid in this region. A naive NLO calculation shown in Ref. [36] predicts that the spectral function in this region is exponentially suppressed. To model this suppression within the Schrödinger formalism itself, we multiply the imaginary part of the potential by \(e^{-|\omega -2\,M|/T}\) for \(\omega <2\,M\). Moreover, at small spatial distances, \(r \ll 1/m_{\textrm{D}}\), the thermal potential is replaced by a vacuum one [36].

Fig. 3
figure 3

The comparison between perturbative and model spectral functions. The dotted lines represent the perturbative spectral functions, while the solid lines illustrate their modifications as per Eq. (12), for both charmonium (on the left) and bottomonium (on the right). The solid curves correspond to the imaginary-time correlators, labeled as “\(x=\hbox {mod}\)” in Fig. 4, which exhibit good agreement with the lattice data for all distances, except at small distances where the lattice cut-off effects dominate due to the finite lattice spacing

Full size image

After computing both thermal and vacuum contributions to the spectral function, we combine them as

$$\begin{aligned} \rho _{\textrm{PS}}^{\textrm{pert}}(\omega )= A^{\textrm{match}}\rho ^{\textrm{pNRQCD}}_{\textrm{PS}} (\omega )\theta (\omega ^{\textrm{match}}-\omega ) +\rho ^{\textrm{vac}}_{\textrm{PS}}(\omega )\theta (\omega -\omega ^{\textrm{match}})\;. \end{aligned}$$
(11)

In Eq. (11) we introduce a multiplicative factor \(A^{\textrm{match}}\) to the thermal part and it is determined such that both thermal and vacuum parts of the spectral function are connected smoothly at some \(\omega =\omega ^{\textrm{match}}\), where \(2M<\omega ^{\textrm{match}}<3M\). The resulting pseudoscalar spectral functions at various temperatures for both charmonium and bottomonium are shown in Fig. 3. Unlike charmonium, the bottomonium has a resonance peak that is visible at temperatures 220 MeV and 251 MeV. The resonance peak for bottomonium hints the existence of bound state at these given temperatures, however the peak might get broaden at higher temperatures.

5 Modelling and Spectral Reconstruction

In this section, we discuss the process of modeling perturbative spectral functions. We introduce two parameters, A and B, in the perturbative spectral function shown in Fig. 3. Parameter A will account for the normalization of the correlator data, while parameter B will allow for the adjustment of thermal mass shift because pole masses are poorly determined in perturbation theory. As a results the model spectral function gets the following form:

$$\begin{aligned} \rho ^{\textrm{mod}}_{\textrm{PS}}(\omega )=A\rho ^{\textrm{pert}}_{\textrm{PS}}(\omega -B)\;. \end{aligned}$$
(12)

Before fitting the lattice data with the above ansatz, it is important to note that the pseudoscalar correlator can be multiplicatively renormalized, which only affects the fit parameter A. Since our analysis involves correlators at different temperatures obtained from a single lattice spacing, the renormalization constant will be the same for all temperatures. Therefore, we do not need to renormalize the lattice correlator for our present study. Nevertheless, to avoid obtaining large values of A after fitting, we normalize the lattice correlator such that it matches the perturbative correlator at a very short distance of \(\tau T = 0.107\) at a temperature \(T=220\) MeV.

We perform uncorrelated \(\chi ^2\) fit of this normalized lattice data with the ansatz in Eq. (12). The short distance part of the lattice correlator are largely affected by lattice artifacts. As a result, we fit the data starting from \(\tau T\)=0.25 for \(T=220\) MeV and \(\tau T\)=0.21 for \(T=251\) MeV. The resulting fit parameters, A and B, for two distinct temperatures, are presented in Table 3. Preliminary results of the model spectral function \(\rho ^{\textrm{mod}}\) and the lattice correlator data, are shown in Fig. 4. The important observation is the shift of the thresholds towards larger energy in the model spectral function. Overall, the observed shift is perfectly admissible and could be due to the uncertainties in the determination of the pole mass or the thermal mass corrections. It is important to take these uncertainties into account when interpreting the results.

Table 3 Estimates of the best fit parameters at two temperatures according to Eq. (12) for charmonium and bottomonium
Full size table
Fig. 4
figure 4

The data points correspond to the lattice correlators, while the lines depict the model correlators obtained from Eq. (1) for charmonium (on the left) and bottomonium (on the right)

Full size image

6 Conclusion and Outlook

In this work we have presented preliminary results of the spectral reconstruction in the pseudoscalar channel for quarkonium, using an ansatz for the fitting procedure based on perturbation theory. We solved the Schrödinger equation for the thermal part of the spectral function around the threshold, using the perturbative finite temperature potential. For energies well above the threshold, we considered the vacuum part and matched both parts smoothly. Our analysis was limited to two temperatures, namely \(T=220\) MeV and \(T=251\) MeV, and we obtained lattice correlators using clover-improved Wilson fermions measured on large \(N_f=2+1\) HISQ gauge field configurations generated at gauge coupling \(\beta =8.249\). we have tuned the quark masses so that the experimental meson spectrum was reproduced on the lattice within errors.

Fitting the spectral function model to the lattice data, we observe a good description of the lattice correlators suggesting that no resonance peaks are needed to describe the charmonium lattice data at these temperatures, while for bottomonium thermally broadened resonance peaks persist at the analyzed temperatures. While our model describes the lattice data well for large time-slice separations, we observed a discrepancy at small distances, which probably is due to the cut-off effects. To estimate these effects, we are in the process of adding more lattice spacings followed by a continuum extrapolation. The non-perturbative effects are important in the temperature range we are considering in this work. Therefore, to improve the model of the spectral function for this range in the future, we will estimate the thermal part of the spectral function using non-perturbatively determined thermal potential from the lattice [51,52,53]. Additionally, we plan to estimate the pole mass, which affects the thermal part of the spectral function and is responsible for the peak position, from the zero-temperature Cornell potential [54]. This will help us better understand the thermal mass shift.

Furthermore, we plan to include lattices at physical quark masses and to consider the vector channel of the quarkonium spectral function, from which we will be able to estimate the heavy quark diffusion coefficients. All data from our calculations, presented in the figures of this paper, can be found in [55].