<figure><img src="image_1.png" /></figure>

# Bottomonium evolution with in-medium heavy quark potential from lattice QCD

Ge Chen, Baoyi Chen, Jiaxing Zhao

► To cite this version:

Ge Chen, Baoyi Chen, Jiaxing Zhao. Bottomonium evolution with in-medium heavy quark poten-
tial from lattice QCD. European Physical Journal C: Particles and Fields, 2024, 84 (8), pp.869.
10.1140/epjc/s10052-024-13266-8 . hal-04483432

HAL Id: hal-04483432

https://hal.science/hal-04483432v1

Submitted on 15 Jan 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Bottomonium evolution with in-medium heavy quark potential from lattice QCD

Ge Chen<fnref n="1" />, Baoyi Chen<fnref n="2" />, Jiaxing Zhao<fnref n="3" />,ᵃ ID

Received: 2 April 2024 / Accepted: 20 August 2024 / Published online: 30 August 2024
© The Author(s) 2024

**Abstract** The static properties and dynamic evolution of bottomonium states in a hot QCD medium are investigated through the Schrödinger equation with complex heavy quark potentials, which have been presented recently in a lattice QCD study and with three different extractions. This approach builds a direct connection between the in-medium heavy quark potentials from the lattice QCD to the experimental observables. The yields and nuclear modification factors $R_{AA}$ of bottomonium in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV are calculated in this work. Our results show a large suppression of the bottomonium yield in heavy ion collisions due to the large imaginary potential. To understand bottomonium $R_{AA}$ based on lattice QCD potentials, we propose a formation time for bottomonium states and find that experimental data can be well explained with the heavy quark potential extracted by the Padé fit, which shows no color screening in the real part potential.

It is widely believed that relativistic heavy-ion collisions, performed at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC), can generate the deconfined Quantum Chromodynamic (QCD) matter called Quark-Gluon Plasma (QGP) [1,2]. Extensive research has been conducted to study the properties of the QGP at different temperatures and baryon chemical potentials. Quarkonia, the bound state of heavy quark and its antiquark, is an excellent probe to study the QGP properties. Besides, quarkonia can be used to examine the QCD at finite temperatures. Due to the hierarchy of energy scales $m_Q \gg m_{Qv} \gg m_{Qv}^2$ of heavy quarks in vacuum, the quarkonia can be described by, for example, the non-relativistic QCD (NRQCD) and potential non-relativistic QCD (pNRQCD) theories [3–5]. The equation of motion in pNRQCD returns to a Schrödinger-like equation with the transition between

color singlet state and octet state [5], where the potential at the leading order for singlet state reads $V(r) = -C_F \alpha_s / r$ with $C_F = (N_c^2 - 1) / (2N_c)$ and the number of colors $N_c = 3$. Furthermore, the non-perturbative lattice simulations show a linear potential between a heavy quark-antiquark pair at long range, and the full range potential is very close to the Cornell potential [6]. Employing the two-body Schrödinger equation with the Cornell potential, the quarkonium mass spectra can be described successfully, see the review paper [7].

In the hot medium, the perturbative investigation of quarkonium properties depends on the relationship of new scales, e.g. temperature $T$, Debye screening mass $m_D$. The previous perturbative studies indicate that in most cases, quarkonium properties in the hot medium are determined by a finite-temperature potential, which has an imaginary part [8,9]. The real part of the potential reveals the color screening effect in the hot medium [10], while the imaginary part is introduced by the Landau damping or color singlet to octet transition [8]. At extremely high temperatures, the Landau damping becomes dominant, and the potential can be obtained by the Hard-Thermal Loop (HTL) resummed perturbation [11,12]. In the strong coupling regime, various methods such as lattice QCD [13–17] and the T-matrix method [18] have been employed to investigate the heavy quark potential.

With quarkonium thermal properties, their production in heavy ion collisions has been approached through various methods. The classical transport approach [18–26], quantum method such as time-dependent Schrödinger equation [27–31], the improved Remler equation [32,33], and open quantum system approach [34–43]. The final yield and transverse momentum spectrum reveal the in-medium binding energy and thermal width of quarkonium, which are encoded in the real and imaginary part of the potential. In the framework of the time-dependent Schrödinger equation, the color screen-

<footnote n="1">¹ Department of Physics and Astronomy, UCL, Gower St, London WC1E 6BT, England</footnote>

<footnote n="2">² Department of Physics, Tianjin University, Tianjin 300354, China</footnote>

<footnote n="3">³ SUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS, 4 rue Alfred Kastler, 44307 Nantes cedex 3, France</footnote>

<footnote>ᵃ e-mail: zhao-jx15@tsinghua.org.cn (corresponding author)</footnote>
ing effect encoded in the real part of the potential results in transitions between different bound states or scattering states, while the imaginary part of the potential suppresses the wave function of singlet states which can be attributed to the transition from singlet to octet states. As components of octet states are usually not included in the Schrödinger equation, those singlet-octet state transitions reduce the normalization of the wave function. In previous studies, the real part of the in-medium potential is suggested to be located between the free energy and the internal energy, both of which show a clear color screening effect as temperature increases, while the imaginary potential is given by the HTL or quenched lattice QCD [28–31]. However, recent lattice QCD calculations with dynamic quarks suggest that there are no color screening effects for the real part heavy quark potential, even at high temperatures, e.g. $T \le 3T_c$ with the critical temperature $T_c$ and $r \le 1$ fm/c [17,44], which is larger than bottomonium averaged radius [7]. In the meantime, the imaginary part is much larger than that from the HTL and quenched lattice QCD. This unscreened or less screened potential has been shown in other studies. The heavy quark potential has been extracted via a machine learning method from the bottomonium masses and widths, which are given by the lattice NRQCD [45]. The result shows the color screening is negligible. Compared to the perturbative results (such as HTL), it seems that this non-screening potential comes from the non-perturbative contribution. In this sense, the Gribov-Zwanziger approach, which gives an infrared improved gluon propagator, has been used to study the heavy quark potential [46,47]. The result shows that the in-medium heavy quark potential has a less color screening and a large imaginary part.

In this paper, we utilize the heavy quark finite-temperature potential, which has been obtained recently in lattice QCD, to study the bottomonium static properties and dynamical evolution in heavy ion collisions.

In lattice QCD, the heavy quark potential can be extracted from the spectrum function. To constrain the spectral function from limited data on Euclidean time correlation functions, one needs to assume some functional form for it. The spectral function has a low-frequency tail, a temperature-independent high-frequency term, and a peak structure at intermediate frequency regions. Recent studies employed several different functional forms to extract the spectrum function, such as the Gaussian fit (a Gaussian shape for the peak), the HTL fit (the peak position and peak width are considered as the real and imaginary part of thermal static energy, which is given by HTL method), and Padé fit (a model-independent approach and based on the Padé rational approximation). The results indicate that the extraction method significantly influences the resulting heavy quark potential [17]. In the following, we will introduce these three different potentials obtained by the above-mentioned extraction methods. However, the lattice results are discrete points in coordinate space $r$ and

with several given temperatures, while we need a temperature and coordinate $r$ smoothly dependent potential for the time-dependent Schrödinger equation. In this sense, we fitted the lattice data with three different pasteurization forms. The first one is extracted by the Gaussian fits method. The complex potential can be separated as a real part and an imaginary part, $V(r, T) = V_R(r, T) - iV_I(r, T)$. The real and imaginary parts of the potential are shown on the left of Fig. 1. The real part potential shows non-screening and the temperature-reduced imaginary part potential has the same trend. We parametrize the potential with the form,

$$ 
\begin{aligned}
V_R^{\text{Gauss}}(r, T) &= -\frac{\alpha_1}{r} + \sigma_1 r, \\
V_I^{\text{Gauss}}(r, T) &= T[(rT)^{\beta_1} + \gamma_1(rT)],
\end{aligned} 
\quad (1) $$

where the parameters $\alpha_1 = 0.3805$ and $\sigma_1 = 0.22 \text{ GeV}^2$, $\beta_1 = 1.2$ and $\gamma_1 = 0.54$. Besides, the real and imaginary potential obtained via the HTL extraction [17] show a huge difference where an obvious screening effect appears. The lattice data and the fitting curves are shown in the middle of Fig. 1. The parameterization form can be expressed as,

$$ 
\begin{aligned}
V_R^{\text{HTL}}(r, T) &= -\alpha_2 \left( m_{\text{DR}}^{\text{HTL}} + \frac{e^{-m_{\text{DR}}^{\text{HTL}}r}}{r} \right) \\
&\quad + \frac{2\sigma_2}{m_{\text{DR}}^{\text{HTL}}} - \frac{e^{-m_{\text{DR}}^{\text{HTL}}r}(2 + m_{\text{DR}}^{\text{HTL}}r)\sigma_2}{m_{\text{DR}}^{\text{HTL}}}, \\
V_I^{\text{HTL}}(r, T) &= \gamma_2\alpha_2 T \phi(m_{\text{DI}}^{\text{HTL}}r),
\end{aligned} 
\quad (2) $$

where $\alpha_2 = 0.3805$, $\sigma_2 = 0.22 \text{ GeV}^2$, $\gamma_2 = 100$. The Debye screening mass in real ($m_{\text{DR}}^{\text{HTL}}$) and imaginary part ($m_{\text{DI}}^{\text{HTL}}$) can be obtained by fitting the lattice data. The definition of $\phi$ is

$$ \phi(x) = 2 \int_{0}^{\infty} dz \frac{z}{(z^2 + 1)^2} \left( 1 - \frac{\sin(xz)}{xz} \right). \quad (3) $$

The third one is obtained via the Padé fit as shown in the right of Fig. 1, which also has a non-screening real part potential. The parameterization of $V(r, T)$ can be expressed as,

$$ 
\begin{aligned}
V_R^{\text{Pade}}(r, T) &= -\frac{\alpha_3}{r} + \sigma_3 r, \\
V_I^{\text{Pade}}(r, T) &= \gamma_3\alpha_3 T \phi(m_{\text{DI}}^{\text{Pade}}r),
\end{aligned} 
\quad (4) $$

where $\alpha_3 = 0.41$, $\sigma_3 = 0.24 \text{ GeV}^2$, and $\gamma_3 = 60$. The Debye screening mass in the imaginary part $m_{\text{DI}}^{\text{Pade}}$ can be obtained from the lattice data.

The static properties of the bottomonium states in a hot medium can be calculated by the Schrödinger equation with the above-mentioned finite-temperature potential at a given temperature. Because the potential is angle-independent, the center-of-mass motion can be separated and the reduced two-
<figure><img src="image_3.png" /><figcaption>Fig. 1 The real and imaginary part potential between heavy quark Q and $\bar{Q}$ at finite temperature. From left to right are the potential obtained from the Gaussian fits, the HTL fits, and the Padé fits. The data are from the lattice QCD simulation [17]. The curves are the fitting results with Eq. (1), (2), and (4)</figcaption></figure>

body Schrödinger equation reads,

$$ \left[-\frac{\nabla^2}{2m_\mu} + V_R(r, T) - iV_I(r, T)\right] \Psi(r) = E_r \Psi(r), \quad (5) $$

where $m_\mu = (m_b m_{\bar{b}})/(m_b + m_{\bar{b}}) = m_b/2$ is the reduced mass and bottom quark mass $m_b = 4.7$ GeV, $\psi$ is relative wavefunction of bottomonium states and $E_r$ is the eigen energies. With the introduction of this complex potential, both the energy eigenvalues and wavefunctions acquire complex values. The real part of the energy eigenvalue determines the mass of the bottomonium state: $M(T) = 2m_b + \text{Re}[E_r(T)]$, while the imaginary part corresponds to the thermal width $\Gamma(T) = -\text{Im}[E_r(T)]$. The averaged radius of the bottomonium state is defined as $\langle r \rangle = \int \Psi\Psi^* |r| d^3r$. Figure 2 illustrates the masses, thermal widths, and radii of the bottomonium states up to the 3S state. It is observed that the potentials obtained from the Gaussian fits and Padé fits, which lack the screening effect, have almost no impact on the real part of the potential. Consequently, the masses using these two potential models remain nearly constant and equivalent to the vacuum values. Conversely, the HTL fits yield a potential with strong screening, resulting in a decrease in mass until the bound

state disappears and an increase in the averaged radius. As for the thermal widths, there are noticeable differences among the three potential models. The larger imaginary potential derived from the Gaussian fits gives rise to a wider thermal width compared to the other two models.

With these static properties in the hot medium, we now come to study the dynamic evolution and production of bottomonium states in relativistic heavy ion collisions. The quark-gluon plasma has been observed to be a strongly-coupled medium created in relativistic heavy ion collisions. Its dynamical expansion can be effectively simulated using hydrodynamic equations in conjunction with an appropriate equation of state (EoS). In this study, we utilized the MUSIC package [48,49] to simulate the expansion of hot medium and adopted the parametrization of the lattice QCD EoS with a smooth crossover between the QGP and the hadron resonance gas around critical temperature $T_c = 170$ MeV [50,51]. The initial temperature profiles of the hot medium were determined by analyzing the final multiplicity of charged hadrons. For the most central collisions and the central rapidity region, the maximum initial temperature of the medium was calculated to be $T_0(\tau_0 = 0.6, \mathbf{x}_T = 0) = 510$ MeV at the center
<figure><img src="image_4.png" /><figcaption>Fig. 2 The masses, widths, and radii of various bottomonium states under different finite-temperature potentials are depicted in Fig. 1. The potential obtained from the Gaussian fits is represented by the long-dashed lines, the HTL fits by the short-dashed lines, and the Padé fits by the solid lines</figcaption></figure>

of the medium in $\sqrt{s_{NN}}$ = 5.02 TeV Pb-Pb collisions [52]. Here, $\tau_0$ = 0.6 fm/c refers to the start time of the hydrodynamics. An effective shear viscosity $\eta/s$ = 0.08 and a zero bulk viscosity are chosen in this study [51].

The crucial point is the formation time $\tau_\gamma$, also called the decoherence time of the bottomonium states, which is a quantum effect and still not clear.

The heavy quark pair carries a narrow wave packet after it is produced at $\tau_{b\bar{b}}$ ≈ 1/(2m_b), and this wave packet evolves into a broader wave package. We neglect the detailed evolution of the wave packet in the QCD medium at pre-equilibrium stage and assume that it has evolved to a certain bottomonium eigenstate at the time $\tau_\gamma$. Different bottomonium states may have different formation times $\tau_\gamma$ that can be estimated via $\tau_\gamma \propto 1/E_B$ with the binding energy of the bottomonium state. So, a longer formation time is expected for higher excited states. We employ the same value of the formation time for different eigenstates for simplicity. There are no bottomonium states before $\tau_\gamma$ instead of a wave packet which is a superposition state of all bottomonium states. The Schrödinger equation starts at $\tau > \tau_\gamma$, where the initial conditions are taken as different bottomonium eigenstates.

Because the mass of bottom quarks is very large, the relative motion of bottom quarks in the bottomonium can be treated in a nonrelativistic way. One can therefore employ the time-dependent Schrödinger equation to describe the evolutions of bottomonium wave functions with in-medium potentials [28, 31, 53, 54]. The in-medium potential is almost isotropic when neglecting the shear-viscosity correction. So, the center-of-mass momentum is conserved and the center-of-mass motion can be separated. Furthermore, the wave-function of the relative motion can be factorized into a radial part and angular part, $\Psi(r, t) = R_{nl}(r, t)Y_{lm}(\theta, \phi)$. The radial wavefunction of state $(n, l)$ satisfies the Schrödinger equation,

$$i \frac{\partial}{\partial t} \psi_{nl}(r, t) = \hat{H} \psi_{nl}(r, t), \\ \hat{H} = -\frac{1}{2m_{\mu}} \frac{\partial^2}{\partial r^2} + V_R(r, T) - iV_I(r, T) + \frac{L(L+1)}{2m_{\mu}r^2}, \quad (6)$$

where $\psi_{nl}(r, t)$ is defined as $\psi_{nl}(r, t) = r R_{nl}(r, t)$. L is the quantum number of the angular momentum. $m_\mu$ is the reduced mass. The in-medium potential can be described by the three aforementioned parametrizations. As a result of the medium-modified potentials, there are internal evo-
lutions within the $b\bar{b}$ wave package, leading to transitions between different bottomonium eigenstates. To study these evolutions, we employ the Crank-Nicolson method to numerically solve the time-dependent Schrödinger equation.

Because the total momentum is conserved, $b\bar{b}$ wave package performs a uniform linear motion with the given initial momentum in the hot medium. In heavy ion collisions, there is more than one $b\bar{b}$ package primordially produced at different positions with different total momentum. $b\bar{b}$ packages in nuclear collisions can be treated as a superposition of the distribution in proton-proton (p-p) collisions [30]. The initial spatial distribution of the $b\bar{b}$ wave package is proportional to the density of binary collisions. The momentum distribution of the total momentum can be estimated by the observed momentum spectrum of $\Upsilon$ in p-p collisions and can be fitted via the formula,

$$ f_{pp}^{b\bar{b}}(\mathbf{p}_T) = \frac{d\tilde{N}_{pp}^{b\bar{b}}}{2\pi p_T dp_T} = \frac{(n-1)}{\pi(n-2)\langle p_T^2 \rangle_{pp}} \left[ 1 + \frac{p_T^2}{(n-2)\langle p_T^2 \rangle_{pp}} \right]^{-n} \quad (7) $$

with $n = 2.5$ and the mean transverse momentum square at the central rapidity is fitted to be $\langle p_T^2 \rangle_{pp} = 80$ (GeV/c)$^2$ at 5.02 TeV [30] based on the measurements of $\Upsilon(1S)$ in p-p collisions [57,58].

When the $b\bar{b}$ wave package moves out of the QGP medium, the final fraction of bottomonium eigenstates labeled with the radial and angular quantum number $(n, l)$ is calculated by projecting the wave package to the corresponding vacuum eigenstate wave function $|c_{nl}|^2 = |\langle\psi(r, t)|\psi_{nl}^{\text{vacuum}}(r)\rangle|^2$. The yield of bottomonium eigenstate $(n, l)$ in nuclear collisions can be obtained by averaging over the spatial and momentum space [30].

$$ \langle |c_{nl}(t)|^2 \rangle = \frac{\int d\mathbf{x}_T d\mathbf{p}_T |c_{nl}(t, \mathbf{x}_T, \mathbf{p}_T)|^2 f_{pp}^{b\bar{b}}(\mathbf{p}_T, \mathbf{b}) n_{\text{coll}}(\mathbf{x}_T) \mathcal{R}(\mathbf{x}_T)}{\int d\mathbf{x}_T d\mathbf{p}_T f_{pp}^{b\bar{b}}(\mathbf{p}_T) n_{\text{coll}}(\mathbf{x}_T, \mathbf{b})} \quad (8) $$

where $n_{\text{coll}}(\mathbf{x}_T, \mathbf{b}) = \sigma_{\text{inel}} T_A(\mathbf{x}_T + \mathbf{b}/2) T_A(\mathbf{x}_T - \mathbf{b}/2)$ is the density of binary collisions with impact parameter $\mathbf{b}$. $T_{A(B)}$ is the thickness functions [59] of nucleus $A(B)$ and $\sigma_{\text{inel}}$ is the inelastic scattering cross-section between nucleons. For Pb-Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV, $\sigma_{\text{inel}} = 68$ mb [60]. $\mathcal{R}$ is the shadowing effect, which is calculated from the EPS09 NLO package [61],

The number of initially produced bottomonium states will be changed in the hot medium due to either the transition to other states or the dissociation, which describes the transition from the color singlet state to the octet state and is encoded in the imaginary potential. The nuclear modification factor

reveals the number change directly. Due to the decays from various states, the final nuclear modification factor for $\Upsilon(nS)$ can be expressed as follows,

$$ R_{AA}(\Upsilon(ns)) = \frac{\sum_{nl} \langle |c_{nl}(t)|^2 \rangle \sigma_{pp}^{nl} \mathcal{B}_{nl \to ns}}{\sum_{nl} \langle |c_{nl}(t_0)|^2 \rangle \sigma_{pp}^{nl} \mathcal{B}_{nl \to ns}}. \quad (9) $$

where $\langle |c_{nl}(t_0)|^2 \rangle$ in the denominator characterizes the initial production of bottomonium eigenstates without both hot and cold nuclear matter effects, so it can also be calculated with Eq. (8) by substituting $\mathcal{R} = 1$ into the equation. The direct production cross section of different bottomonium eigenstates $\sigma_{pp}^{nl}$ is extracted from the experimental data in p-p collisions, $\sigma_{pp}(1s, 1p, 2s, 2p, 3s) = (37.97, 44.20, 18.27, 37.68, 8.21)$ nb [28,30] without feed-down contribution. While the branching ratio of the decay $\mathcal{B}_{nl \to 1s}$ is taken from Particle Data Group [62]. The decay values from higher to lower states are taken as $\mathcal{B}_{(1p,2s,2p,3s) \to 1s} = (0.247, 0.264, 0.088, 0.065)$, $\mathcal{B}_{(2p,3s) \to 2s} = (0.124, 0.106)$ and $\mathcal{B}_{(2p,3s) \to 1p} = (0.0065, 0)$.

Similar equation is employed to get the prompt $R_{AA}(2S)$ and $R_{AA}(3S)$.

Now let's see the results obtained by taking the three different in-medium potentials as discussed in the last section. First, we assume that all bottomonium eigenstates have been formed at $\tau_\Upsilon = \tau_0 = 0.6$ fm/c, which is the start time of the hydrodynamics, and subsequently, their internal evolution is governed by the time-dependent Schrödinger equation (6) with complex potentials. The nuclear modification factors of bottomonium as a function of the number of participants ($N_{part}$) and transverse momentum ($p_T$) are calculated with three different potentials shown in Fig. 3. Due to the very large imaginary potentials in the Gaussian fit, bottomonium suppression is significant and far below the experimental data points. When the imaginary potential is taken as HTL or Padé fit, the bottomonium suppression becomes a bit weaker, but still underestimates the experimental data, see Fig. 3. With this formation time, the bottomonium stays the longest time in the hot medium and suffers the largest dissociation in the hot medium.

In another limit, all bottomonium states are formed with a large value of formation time $\tau_\Upsilon = 3$ fm/c. Other setups are the same as before. In this case, the dissociation rates of bottomonium states are zero before 3 fm/c even if they are in the QGP medium. The results are shown in Fig. 4. The hot medium suppression becomes weaker for all potentials. We find the calculations with Padé fit can well explain the data points.

In our studies, the real part of the potential does not change the nuclear modification factors of bottomonium. The imaginary part of the potential, which contributes an exponential damping factor in the wave function of $b\bar{b}$, can significantly affect the final results. In our recent framework, the dissocia-
<figure><img src="image_5.png" /><figcaption>**Fig. 3** The nuclear modification factors of bottomonium $Υ(1s, 2s, 3s)$ as a function of the number of participants $N_{part}$ (upper) and the transverse momentum $p_T$ (lower) in the central rapidity of Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The complex potential is taken as The in-medium</figcaption></figure>

potential is taken from the Gaussian fit (left), HTL fit (middle), and the Padé fit (right), respectively. The experimental data are cited from CMS [55,56]

tion from the color singlet state to the octet state is involved in the imaginary potential in the Hamiltonian, while the inverse process from the octet to singlet state has been temporarily neglected, which can be interpreted as regeneration from a correlated $b\bar{b}$ dipole. Now we want to discuss the importance of transition from octet to singlet state. If one starts to evolve the $b\bar{b}$ states from a singlet state, then of course the transition from singlet to octet is dominated while the inverse process is largely suppressed as used in previous studies [42,43]. After a long time evolution, the transition will approach detail balance. However, if one starts to evolve the $b\bar{b}$ states from an octet state, the transition from the octet to the singlet state is definitely important and cannot be neglected.

In summary, we employ the Schrödinger equation with complex potentials obtained recently from lattice QCD to study the static properties and dynamical evolution of bottomonium states in the quark-gluon plasma, created in relativistic heavy ion collisions. The mass and thermal width of

the bottomonium state reveal the real and imaginary parts of the in-medium potential, respectively. The evolution of bottomonium wave functions is studied under different potentials, where the corresponding nuclear modification factors are calculated. Due to the large imaginary potential in Gaussian fit, the nuclear modification factors are reduced to almost zero. While with the other two complex potentials from the HTL fit and the Padé fit, the medium suppression becomes weaker than the case of the Gaussian fit, but still lower than the experimental data. We find with a large formation time, the $R_{AA}$ of bottomonium can be well described with the heavy quark potential extracted by the Padé fit, which shows no color screening in the real part potential. This indicates the importance of the bottomonium formation time. This framework helps to establish the correspondence between the in-medium heavy quark potential from first-principle-based lattice QCD and experimental observables about heavy quarkonium, aiding in the understanding of heavy quark physics.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Same as Fig. 3 but with the formation time of bottomonium $\tau_{\gamma} = 3$ fm/c</figcaption></figure>

**Acknowledgements** This work is supported by the National Natural Science Foundation of China (NSFC) under Grants no. 12175165. JX is supported by the European Union's Horizon 2020 research and innovation program under Grant agreement No. 824093 (STRONG-2020).

**Data Availability Statement** Data will be made available on reasonable request. [Authors' comment: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.]

**Code Availability Statement** Code/software will be made available on reasonable request. [Authors' comment: The code/software generated during and/or analysed during the current study is available from the corresponding author on reasonable request.]

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP³.

## References

1. A. Bazavov et al., Phys. Rev. D **85**, 054503 (2012). https://doi.org/10.1103/PhysRevD.85.054503. arXiv:1111.1710 [hep-lat]

2. W. Busza, K. Rajagopal, W. van der Schee, Annu. Rev. Nucl. Part. Sci. **68**, 339 (2018). https://doi.org/10.1146/annurev-nucl-101917-020852. arXiv:1802.04801 [hep-ph]

3. W.E. Caswell, G.P. Lepage, Phys. Lett. B **167**, 437 (1986). https://doi.org/10.1016/0370-2693(86)91297-9

4. N. Brambilla, A. Pineda, J. Soto, A. Vairo, Nucl. Phys. B **566**, 275 (2000). https://doi.org/10.1016/S0550-3213(99)00693-8. arXiv:hep-ph/9907240

5. N. Brambilla, A. Pineda, J. Soto, A. Vairo, Rev. Mod. Phys. **77**, 1423 (2005). https://doi.org/10.1103/RevModPhys.77.1423. arXiv:hep-ph/0410047

6. T. Kawanai, S. Sasaki, Phys. Rev. D **85**, 091503 (2012). https://doi.org/10.1103/PhysRevD.85.091503. arXiv:1110.0888 [hep-lat]

7. J. Zhao, K. Zhou, S. Chen, P. Zhuang, Prog. Part. Nucl. Phys. **114**, 103801 (2020). https://doi.org/10.1016/j.ppnp.2020.103801. arXiv:2005.08277 [nucl-th]

8. N. Brambilla, J. Ghiglieri, A. Vairo, P. Petreczky, Phys. Rev. D **78**, 014017 (2008). https://doi.org/10.1103/PhysRevD.78.014017. arXiv:0804.0993 [hep-ph]

9. N. Brambilla, M.A. Escobedo, J. Ghiglieri, J. Soto, A. Vairo, JHEP **09**, 038. https://doi.org/10.1007/JHEP09(2010)038. arXiv:1007.4156 [hep-ph]

10. F. Karsch, M.T. Mehr, H. Satz, Z. Phys. C **37**, 617 (1988). https://doi.org/10.1007/BF01549722
11. M. Laine, O. Philipsen, P. Romatschke, M. Tassler, JHEP 03, 054. https://doi.org/10.1088/1126-6708/2007/03/054. arXiv:hep-ph/0611300

12. A. Beraudo, J.P. Blaizot, C. Ratti, Nucl. Phys. A 806, 312 (2008). https://doi.org/10.1016/j.nuclphysa.2008.03.001. arXiv:0712.4394 [nucl-th]

13. A. Rothkopf, T. Hatsuda, S. Sasaki, Phys. Rev. Lett. 108, 162001 (2012). https://doi.org/10.1103/PhysRevLett.108.162001. arXiv:1108.1579 [hep-lat]

14. D. Lafferty, A. Rothkopf, Phys. Rev. D 101, 056010 (2020). https://doi.org/10.1103/PhysRevD.101.056010. arXiv:1906.00035 [hep-ph]

15. A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo, J.H. Weber (TUMQCD), Phys. Rev. D 98, 054511 (2018). https://doi.org/10.1103/PhysRevD.98.054511. arXiv:1804.10600 [hep-lat]

16. P. Petreczky, J. Weber (TUMQCD), Nucl. Phys. A 967, 592 (2017). https://doi.org/10.1016/j.nuclphysa.2017.04.011. arXiv:1704.08573 [hep-lat]

17. D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf, J.H. Weber (HotQCD), Phys. Rev. D 105, 054513 (2022). https://doi.org/10.1103/PhysRevD.105.054513. arXiv:2110.11659 [hep-lat]

18. S.Y.F. Liu, R. Rapp, Phys. Rev. C 97, 034918 (2018). https://doi.org/10.1103/PhysRevC.97.034918. arXiv:1711.03282 [nucl-th]

19. X. Zhao, R. Rapp, Phys. Rev. C 82, 064905 (2010). https://doi.org/10.1103/PhysRevC.82.064905. arXiv:1008.5328 [hep-ph]

20. M. He, B. Wu, R. Rapp, Phys. Rev. Lett. 128, 162301 (2022). https://doi.org/10.1103/PhysRevLett.128.162301. arXiv:2111.13528 [nucl-th]

21. X. Du, R. Rapp, Phys. Lett. B 834, 137414 (2022). https://doi.org/10.1016/j.physletb.2022.137414. arXiv:2207.00065 [nucl-th]

22. L. Yan, P. Zhuang, N. Xu, Phys. Rev. Lett. 97, 232301 (2006). https://doi.org/10.1103/PhysRevLett.97.232301. arXiv:nucl-th/0608010

23. K. Zhou, N. Xu, Z. Xu, P. Zhuang, Phys. Rev. C 89, 054911 (2014). https://doi.org/10.1103/PhysRevC.89.054911. arXiv:1401.5845 [nucl-th]

24. B. Chen, Chin. Phys. C 43, 124101 (2019). https://doi.org/10.1088/1674-1137/43/12/124101. arXiv:1811.11393 [nucl-th]

25. X. Yao, T. Mehen, Phys. Rev. D 99, 096028 (2019). https://doi.org/10.1103/PhysRevD.99.096028. arXiv:1811.07027 [hep-ph]

26. X. Yao, B. Müller, Phys. Rev. D 100, 014008 (2019). https://doi.org/10.1103/PhysRevD.100.014008. arXiv:1811.09644 [hep-ph]

27. R. Katz, P.B. Gossiaux, Ann. Phys. 368, 267 (2016). https://doi.org/10.1016/j.aop.2016.02.005. arXiv:1504.08087 [quant-ph]

28. A. Islam, M. Strickland, JHEP 21, 235. https://doi.org/10.1007/JHEP03(2021)235. arXiv:2010.05457 [hep-ph]

29. A. Islam, M. Strickland, Phys. Lett. B 811, 135949 (2020). https://doi.org/10.1016/j.physletb.2020.135949. arXiv:2007.10211 [hep-ph]

30. L. Wen, B. Chen, Phys. Lett. B 839, 137774 (2023). https://doi.org/10.1016/j.physletb.2023.137774. arXiv:2208.10050 [nucl-th]

31. L. Wen, X. Du, S. Shi, B. Chen, Chin. Phys. C 46, 114102 (2022). https://doi.org/10.1088/1674-1137/ac7fe6. arXiv:2205.07520 [nucl-th]

32. D.Y.A. Villar, J. Zhao, J. Aichelin, P.B. Gossiaux, Phys. Rev. C 107, 054913 (2023). https://doi.org/10.1103/PhysRevC.107.054913. arXiv:2206.01308 [nucl-th]

33. T. Song, J. Aichelin, J. Zhao, P.B. Gossiaux, E. Bratkovskaya, Phys. Rev. C 108, 054908 (2023). https://doi.org/10.1103/PhysRevC.108.054908. arXiv:2305.10750 [nucl-th]

34. N. Borghini, C. Gombeaud, Eur. Phys. J. C 72, 2000 (2012). https://doi.org/10.1140/epjc/s10052-012-2000-7. arXiv:1109.4271 [nucl-th]

35. J.-P. Blaizot, D. De Boni, P. Faccioli, G. Garberoglio, Nucl. Phys. A 946, 49 (2016). https://doi.org/10.1016/j.nuclphysa.2015.10.011. arXiv:1503.03857 [nucl-th]

36. J.-P. Blaizot, M.A. Escobedo, JHEP 06, 034. https://doi.org/10.1007/JHEP06(2018)034. arXiv:1711.10812 [hep-ph]

37. J.-P. Blaizot, M.A. Escobedo, Phys. Rev. D 98, 074007 (2018). https://doi.org/10.1103/PhysRevD.98.074007. arXiv:1803.07996 [hep-ph]

38. Y. Akamatsu, A. Rothkopf, Phys. Rev. D 85, 105011 (2012). https://doi.org/10.1103/PhysRevD.85.105011. arXiv:1110.1203 [hep-ph]

39. Y. Akamatsu, Phys. Rev. D 87, 045016 (2013). https://doi.org/10.1103/PhysRevD.87.045016. arXiv:1209.5068 [hep-ph]

40. S. Kajimoto, Y. Akamatsu, M. Asakawa, A. Rothkopf, Phys. Rev. D 97, 014003 (2018). https://doi.org/10.1103/PhysRevD.97.014003. arXiv:1705.03365 [nucl-th]

41. S. Delorme, T. Gousset, R. Katz, P.-B. Gossiaux, EPJ Web Conf. 258, 05009 (2022). https://doi.org/10.1051/epjconf/202225805009

42. N. Brambilla, M.A. Escobedo, J. Soto, A. Vairo, Phys. Rev. D 96, 034021 (2017). https://doi.org/10.1103/PhysRevD.96.034021. arXiv:1612.07248 [hep-ph]

43. N. Brambilla, M.A. Escobedo, A. Islam, M. Strickland, A. Tiwari, A. Vairo, P. Vander Griend, Phys. Rev. D 108, L011502 (2023). https://doi.org/10.1103/PhysRevD.108. L011502. arXiv:2302.11826 [hep-ph]

44. A. Bazavov, D. Hoying, O. Kaczmarek, R.N. Larsen, S. Mukherjee, P. Petreczky, A. Rothkopf, J.H. Weber. (2023). arXiv:2308.16587 [hep-lat]

45. S. Shi, K. Zhou, J. Zhao, S. Mukherjee, P. Zhuang, Phys. Rev. D 105, 014017 (2022). https://doi.org/10.1103/PhysRevD.105.014017. arXiv:2105.07862 [hep-ph]

46. W. Wu, G. Huang, J. Zhao, P. Zhuang, Phys. Rev. D 107, 114033 (2023). https://doi.org/10.1103/PhysRevD.107.114033. arXiv:2207.05402 [hep-ph]

47. M. Debnath, R. Ghosh, N. Haque. (2023). arXiv:2305.16250 [hep-ph]

48. B. Schenke, S. Jeon, C. Gale, Phys. Rev. Lett. 106, 042301 (2011). https://doi.org/10.1103/PhysRevLett.106.042301. arXiv:1009.3244 [hep-ph]

49. B. Schenke, S. Jeon, C. Gale, Phys. Rev. C 82, 014903 (2010). https://doi.org/10.1103/PhysRevC.82.014903. arXiv:1004.1408 [hep-ph]

50. A. Bazavov et al. (HotQCD), Phys. Rev. D 90, 094503 (2014). https://doi.org/10.1103/PhysRevD.90.094503. arXiv:1407.6387 [hep-lat]

51. J.E. Bernhard, J.S. Moreland, S.A. Bass, J. Liu, U. Heinz, Phys. Rev. C 94, 024907 (2016). https://doi.org/10.1103/PhysRevC.94.024907. arXiv:1605.03954 [nucl-th]

52. W. Zhao, H.-J. Xu, H. Song, Eur. Phys. J. C 77, 645 (2017). https://doi.org/10.1140/epjc/s10052-017-5186-x. arXiv:1703.10792 [nucl-th]

53. B. Krouppa, R. Ryblewski, M. Strickland, Phys. Rev. C 92, 061901 (2015). https://doi.org/10.1103/PhysRevC.92.061901. arXiv:1507.03951 [hep-ph]

54. N. Brambilla, M.A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend, J.H. Weber, JHEP 05, 136. https://doi.org/10.1007/JHEP05(2021)136. arXiv:2012.01240 [hep-ph]

55. A.M. Sirunyan et al. (CMS), Phys. Lett. B 790, 270 (2019). https://doi.org/10.1016/j.physletb.2019.01.006. arXiv:1805.09215 [hep-ex]

56. A. Tumasyan et al. (CMS). (2023). arXiv:2303.17026 [hep-ex]

57. R. Aaij et al. (LHCb), Eur. Phys. J. C 74, 2835 (2014). https://doi.org/10.1140/epjc/s10052-014-2835-1. arXiv:1402.2539 [hep-ex]

58. R. Aaij et al. (LHCb), Eur. Phys. J. C 72, 2025 (2012). https://doi.org/10.1140/epjc/s10052-012-2025-y. arXiv:1202.6579 [hep-ex]
59. M.L. Miller, K. Reygers, S.J. Sanders, P. Steinberg, Annu. Rev. Nucl. Part. Sci. **57**, 205 (2007). https://doi.org/10.1146/annurev.nucl.57.090506.123020. arXiv:nucl-ex/0701025

60. B. Abelev et al. (ALICE), Eur. Phys. J. C **73**, 2456 (2013). https://doi.org/10.1140/epjc/s10052-013-2456-0. arXiv:1208.4968 [hep-ex]

61. K.J. Eskola, H. Paukkunen, C.A. Salgado, JHEP **04**, 065. https://doi.org/10.1088/1126-6708/2009/04/065. arXiv:0902.4154 [hep-ph]

62. M. Tanabashi et al. (Particle Data Group), Phys. Rev. D **98**, 030001 (2018). https://doi.org/10.1103/PhysRevD.98.030001