Heavy quarkonia in a bulk viscous medium

We study the properties of heavy quarkonia in a quark-gluon plasma in the presence of bulk viscous effects. Within the hard thermal loop approximation at one-loop, the dielectric permittivity of quark-gluon plasma is computed, where the bulk viscous effect enters through the deformation of the distribution functions of thermal quarks and gluons. Based on the modified dielectric permittivity, we compute the in-medium heavy quark potential, that includes non-pertubative string-like terms as well as the perturbative Coulombic term. We discuss how the bulk viscous effect modifies the real and imaginary parts of the in-medium potential. Several prescriptions are examined as to how to include the string-like non-perturbative potentials. Using the deformed potential, we compute the wave functions, binding energies, and decay widths of heavy quarkonia in a bulk viscous medium, and study their sensitivity to the strength of the bulk viscous effect. An estimate of the melting temperatures is given.


Introduction
The relativistic heavy-ion collisions provide us with a unique opportunity to experimentally study the strongly interacting matter in extreme conditions. The currently ongoing experimental programs at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN aim at revealing the properties of the quark gluon plasma (QGP), which is expected to appear at high temperatures. At sufficiently high temperatures, a QGP behaves as a weakly interacting gas of quarks and gluons, which can be understood using hard thermal loop (HTL) resummation [1][2][3][4]. Such a description has been successful in describing the thermodynamics of the QGP even close to the crossover temperature [5][6][7].
Heavy quarkonium states have been a useful probe of the surrounding thermal medium. In the vacuum, they are reliably described in terms of non-relativistic potential models [8,9] using the Cornell potential [10,11]. A QGP medium exhibits the screening of static colorelectric fields and that would result in the melting of heavy quarkonia, which was one of the first proposed signals of the formation of a QGP [12]. The potential models have been applied to the study of quarkonia at finite temperatures, the first of such works is done by Karsch, Mehr, and Satz [13]. The meson current correlators and quarkonium spectral functions have been calculated from potential models [14][15][16][17][18][19][20][21] and are compared to the firstprinciple lattice QCD calculations [22][23][24][25][26][27][28]. The appearance of the imaginary part of the potential due do the Landau damping is discussed [29][30][31][32], which has further stimulated the study of complex heavy quarkonia potential from thermal field theories [33][34][35][36] as well as from the lattice QCD [37][38][39]. See Ref. [40] for a recent review.
The motivation for the current work is to understand how the effect of non-equilibrium nature of the fluid is imprinted on the properties of heavy quarkonia. For example, in the early time of a relativistic heavy-ion collision, the longitudinal expansion is stronger than the radial expansion, which would result in an anisotropy of the distribution functions of medium particles in the momentum space. The effect of such momentum-space anisotropies on quarkonia has been discussed in Refs. [33,35,36,[41][42][43][44][45][46][47]. The presence of magnetic fields [48][49][50][51][52][53][54][55][56] or non-zero fluid velocity [57][58][59][60][61][62][63][64] also works as a source of anisotropy. Among such non-equilibrium situations, the role the bulk viscosity is gaining an increasing attention in relation to the beam energy scan program [65], since the bulk viscous effect is expected to be enhanced as the system approaches a critical point [66][67][68].
The goal of this study is to test the sensitivity of quarkonia to the non-equilibrium nature of the fluid, in particular, the bulk viscous corrections. In Ref [69], the dielectric permittivity is computed in the presence of bulk viscous corrections based on the HTL-resummed gluon propagators. The bulk viscous effect enters through the deformation of the distribution functions of thermal particles. The perturbative HTL gluon propagators only gives rise to the Coulombic potential, but non-perturbative string-like contributions have been observed in lattice QCD studies [70][71][72]. There has been several proposed prescriptions as to how to incorporate non-perturbative contributions in the potential [39,44,44,64,73,74]. Among those is an approach based on the linear response theory: the modified string-like potential is obtained by modifying the linear potential using the HTL permittivity that entails the medium effect 1 . In this work, based on the modified dielectric permittivity in the presence of bulk viscous effect, we derive a complex heavy quark potential in such environments. We examine several prescriptions for the introduction of non-perturbative part. We use the modified potential to solve the Schrödinger equation and compute the deformed wave functions, binding energies, and decay widths of heavy quarkonia. We discuss how those physical properties are affected by the bulk viscous effect.
The rest of the article is structured as follows. In Sec. 2, we derive the dielectric permittivity of a thermal medium in the presence of bulk viscous corrections. In Sec. 3, we calculate the complex heavy quark potential based on the modified dielectric permittivity and discuss its properties. In Sec. 4, we show the effect of the bulk viscous corrections on the binding energies and decay widths of quarkonium states, from which melting temperatures are estimated. Section 5 is devoted to summary and discussions.

Dielectric permittivity of a quark gluon plasma with bulk viscous corrections
In order to compute the in-medium potential, we rely on the linear response theory, in which the in-medium properties are encoded in the dielectric permittivity. We here review the derivation of the dielectric permittivity in the HTL approximation in the presence of the bulk viscous correction. When the system is away from thermal equilibrium, the fluctuation-dissipation theorem is violated, which leads to the existence of two different Debye masses. In the current situation, a modified fluctuation-dissipation theorem is found to hold.
In the computation below, the non-equilibrium effect enters through the modification of the distribution function of thermal quarks and gluons, where f 0 (k) is the equilibrium distribution 2 , and the second term is the non-equilibrium correction. In general, non-equilibrium corrections can be anisotropic. Such an anisotropy may be present at the early stage of heavy-ion collisions, where the longitudinal expansion is substantially stronger than the radial expansion. Certain types of the corrections can be regarded as viscous corrections when the anisotropy is weak. In this study, we discuss the effect of the bulk viscosity, and as f 0 we take the thermal equilibrium one, where the correction δ bulk f (k) is isotropic. The specific form of the correction is given later in Eq. (2.7)

Computation of the dielectric permittivity
Let us here derive the dielectric permittivity in the presence of bulk viscous corrections. For this, we compute the gluon self-energies and propagators in the presence of bulk viscous corrections. In the following, we review how to obtain the modified propagators as done in Ref. [69]. The medium quarks are taken to be massless.

Retarded propagator
First, let us look at the retarded self-energy of gluons. For the computation of the potential, we need the temporal component Π R (P ) ≡ Π 00 R (P ) 3 . In the HTL approximation, the oneloop contribution from N f quarks to Π R (P ) is given by [69,75] Π (q) wherek ≡ k/k and f ± (k) are distribution functions for quarks/antiquarks. As long as the HTL approximation is valid, this expression is true even in non-equilibrium situations. In the thermal equilibrium, the distribution functions are given by where µ is the quark chemical potential. The contribution from the gluon loop has the same structure with the Fermi distribution replaced with the Bose one. Including the contribution from quark and gluon loops, the retarded self-energy in the equilibrium is written as where m D is the Debye mass, withμ ≡ µ/T . Now let us introduce the bulk viscous correction. We shall model the correction with the following form [69], where a and Φ are constants and the +(−) sign is for Bose (Fermi) distribution. Φ is a parameter proportional to the bulk viscous pressure (divided by ideal pressure). There are constraints for the parameter a. We need a condition a > 0 so that there will be no IR divergence coming from the correction in the retarded self-energy of gluons (see Eq. (2. 3) with f replaced by the Bose distribution). Otherwise the dominant contribution does not come from k ∼ T and the HTL approximation becomes invalid. The same condition of the absence of IR divergence for the symmetric self-energy (gluon-loop version of Eq. (2.19)) leads to a stronger bound, a > 1/2, from the O(Φ 2 ) contributions. In addition, in order for the two-loop contribution to be negligible, we also need the condition |Φ| g 2 . The bulk pressure δ bulk p is usually negative, which corresponds to Φ < 0, but the sign can be reversed in the presence of shear-bulk coupling [76]. In the present study, we regard Φ as a parameter of either sign.
In the presence of bulk correction (2.7), the retarded self-energy is modified as Π R = Π eq R + δ bulk Π R . The contribution from the quark loop, δ bulk Π (q) R , is given by The correction does not affect the momentum dependence and just modify the Debye mass. Similarly, the contribution from the gluon loop can be computed. The total retarded selfenergy including bulk correction can be written as [69] Π R (P ) = Π eq R + δ bulk Π (q) where m 2 D,R = m 2 D + δm 2 D,R is the modified Debye mass. The correction is written as Here, the dimensionless quantities c q R (a,μ) and c g R (a) are defined by where we take the Bose and Fermi distributions as f 0 for c g R (a) and c q R (a,μ) respectively, and their explicit forms are c g R (a) = 6 π 2 Γ(a + 2)ζ(a + 1), (2.13) where Li n (z) denotes the polylogarithm function. At the vanishing quark chemical potential µ = 0, the quark part c q R simplify to c q R (a,μ = 0) = 12 π 2 (1 − 2 −a )Γ(a + 2)ζ(a + 1). (2.14) We can compute the retarded propagator from the self-energy. In Coulomb gauge, if the distribution function is isotropic, the temporal component of the resummed propagator 4 ,D R (P ) ≡D 00 R (P ) , is independent of the spatial components of the self-energy and propagators 5 . The Dyson-Schwinger equation reads where D R = 1/(p 2 + i sgn(p 0 ) ) is the bare propagator. Using the Π R (P ) obtained above, the temporal component of the resummed retarded propagator is written as For the computation of the potential, we need the static limit p 0 → 0 of the propagator.
To the first order in p 0 , it can be written as where Θ is the step function. The advanced propagator is given by the complex conjugate of the retarded one. 4 Resummed propagators are indicated by characters with bars. 5 In the Coulomb gauge, the bare and resummed propagators satisfy the condition, When the system is isotropic, like in the current case, we have D 0i R = 0.

Symmetric propagator
The symmetric propagator D S (P ) and the self-energy Π S (P ) can be computed in a similar manner. The quark loop contribution to Π S (P ) is written as Adding the quark-loop and gluon-loop contributions, the equilibrium part of Π S (P ) is written using the Debye mass (2.6) as The total symmetric self-energy with the bulk viscous correction is again represented with the modified Debye mass, where m 2 D,S = m 2 D + δm 2 D,S , and the bulk viscous correction δm 2 D,S is given to the first order in Φ by where the dimensionless functions c q S (a,μ) and c g S (a) are defined as Given the symmetric self-energy, we can compute the symmetric propagator. The temporal component of the resummed symmetric propagator in the Coulomb gauge satisfies the following Dyson-Schwinger equation Using the Dyson-Schwinger equation for the retarded and advanced propagators, The first term on the right hand side is in fact zero (note that it is proportional to p 4 δ(p 2 )). Thus, the resummed symmetric propagator can be written asD S =D R Π SDA and in p 0 → 0 limit it is given byD where we have used Eqs. (2.9) and (2.18).

Dielectric permittivity
The dielectric permittivity, (p), is computed as whereD 11 (P ) is the longitudinal component of the 11-part of the resummed gluon propagator. Noting thatD and using Eq. (2.18) and Eq. (2.28), we obtain the expression for the dielectric permittivity in the presence of bulk viscous correction, (2.31) In the limit of vanishing bulk correction, both of m D,R and m D,S approaches m D and the equilibrium expression is reproduced. The dielectric permittivity (2.31) will be used in computing the in-medium heavy quarkonia potential.

Two Debye masses and a modified fluctuation-dissipation theorem
We have learned that the effects of the bulk viscous correction are incorporated in two different Debye masses, m 2 D,R and m 2 D,S , that characterize the retarded (advanced) propagators and the symmetric propagator. The two Debye masses are functions of the bulk viscous correction parameter Φ. We have obtained the expression of the modified dielectric permittivity (2.31), which is the main result of this section.
Let us show the behaviors of the modified Debye masses. The bulk viscous correction can be written as The expression of c R and c S follows from Eqs. (2.10) and (2.22). In Fig. 1, we show the linear coefficients c R (a,μ) and c S (a,μ) as a function of a. Those coefficients are positive in the region of a considered here. Therefore, the bulk viscous correction effectively increases the Debye mass for Φ > 0.
In the absence of bulk viscous corrections, namely in the thermal equilibrium, the fluctuation-dissipation theorem (FDT) holds, 33) and this ensures that the two masses are equal in the absence of bulk viscous corrections. The FDT is violated in non-equilibrium, and m D,R and m D,S can be different. Thus, the difference quantifies the extent of violation of the FDT. In the current situation, in fact, a modified version of FDT holds [77], or equivalently where we have defined a parameter The modified FDT (2.34) holds in the HTL approximation at one-loop and when the distribution function is spherically symmetric. As can be seen in Fig. 1, the symmetric Debye mass is larger than the retarded one for Φ > 0, so λ > 1 in this case.

In-medium potentials in the presence of bulk viscous corrections
In this section, we study how the heavy quarkonia potential is modified in the presence of bulk viscous correction. A heavy-quark potential can be obtained by the Fourier transform of the static gluon propagator. The propagators in the HTL perturbation theory results in the screened Coulombic potential. This potential will be dominant in the high temperature limit, but it does not account for the non-perturbative string-like part, which is responsible for the confinement. At zero or low temperature, many studies has confirmed that the so-called Cornell potential that consists of a Coulombic and string-like parts explain the properties of heavy quarkonia very well. Modeling of non-perturbative effect is important in understanding the "melting" of quarkonia near the crossover temperature T c .
Given the medium property, how to incorporate it to modify the string-like contribution for both of the real and imaginary parts is not unique. There has been a number of proposals as to how to parametrize the in-medium potentials. Here, we shall discuss the prescriptions discussed in Refs. [64,74] and we extend those formulations to introduce the non-equilibrium bulk viscous corrections. We examine how the real and imaginary part of the in-medium potential is affected by this.

Approach based on the linear response
We here take the approach [36] based on the linear response theory. The properties of a thermal medium is encoded in the dielectric permittivity (p). When the linear approximation is justified, the in-medium potential is related to the potential in the vacuum through the permittivity by This relation is true even in a strongly coupled system, as long as the linear approximation to the potential is good. As a vacuum potential, we employ the Cornell potential. Thus, the in-medium heavy-quark potential in the real space can be written as where the Fourier transform of the Cornell potential V Cornell (p) is given by where α ≡ C F α s with C F = (N 2 c − 1)/2N c and σ is the string tension. The parameter σ is determined to reproduce the vacuum quarkonium property.
Let us first look at the real part of the in-medium heavy quark potential. Using Eqs. (2.31) and (3.3), it is computed as where the first term is Coulombic contribution with the Debye screening, and the second term comes from the string-like part of the Cornell potential. In the absence of bulk viscous corrections, this form of potential is derived in Ref. [36] . In the small distance limit, r → 0, it approaches the Cornell potential, V Cornell (r) = −α/r + σr. For the real part of the potential, the bulk viscous correction enters through the modification of the Debye mass m D → m D,R . In Fig. 2, we plot the real part of the potential for different values of Φ (left) and a (right). The modified Debye mass is an increasing function of both Φ and a, and the potential becomes flattened for larger values of those parameters. The heavy quark potential also acquires an imaginary part at finite temperatures. The imaginary part reads where the first term is from the perturbative HTL contribution, and the second term is the string-like contribution. We have plotted these terms separately in Fig. 3, and the total Im VHTL (Φ=0) imaginary part in Fig. 4. Note that the dimensionless parameter λ is defined in Eq. (2.36). The functions φ n (x) and χ(x) are defined by The function φ 2 (x) is a monotonically increasing function that asymptotes φ 2 (0) = 0 and φ 2 (∞) = 1. χ(x) is also monotonically increasing with χ(0) = 0, but is logarithmically divergent at large x 6 . Let us examine the qualitative features of the imaginary part and its bulk viscous correction. We plot the imaginary part of the potential in Fig. 4. As can be seen in the figure, the bulk viscous correction on Im V is different in small r and large r regions. This can be understood as follows: • One consequence of the bulk viscous effect is a shift of the Debye mass, The Debye mass becomes heavier due to the bulk correction for Φ > 0 and φ( m D,R r) increases with the bulk correction. The coefficient of the Coulombic term is αλ, which is also an increasing function of Φ. So, | Im V HTL | is an increasing function of Φ. The perturbative HTL contribution is dominant in the small r region, since φ 2 ( m D,R r) ∼ r 2 ln r and χ ∼ r 2 to the leading order. See the left panel of Fig. 3. Thus, in this region, | Im V | increases as Φ is increased.
• As shown in the right panel of Fig. 3, In the large r region, string-like part dominates the imaginary part. In this region, | Im V string | is suppressed in the presence of bulk correction with Φ > 0. The string part is proportional to the factor λ/ m 2 D,R , which decreases as a function of Φ. Although χ( m D,R r) is an increasing function of Φ, in total the string part decreases as a function of Φ.
• As a result, for the bulk viscous correction for Φ > 0, we observe the enhancement of | Im V | in the small r region, and suppression at large r, as shown in Fig. 4

Introduction of non-perturbative propagator
Let us discuss the prescription given in Ref. [74]. In order to take into account the stringlike behavior of the potential, they have introduced a non-perturbative contribution to the resummed gluon propagator in addition to the HTL contribution as whereD p R = 1//(p 2 − Π R ) is the perturbative HTL contribution. The non-perturbative contribution to the temporal component of the resummed gluon propagator is modeled by the following form,D np where b = 4 and b = 6 is chosen so that the leading contribution to the imaginary part of the potential in the small r limit behaves as r 4 ln r. The mass scale m G is related to the string tension as m 2 G = 2σ/α by matching the small r behavior of the real part of the potential with the Cornell potential. To get the symmetric propagator, the authors used the relationD Using the equilibrium self-energies (2.5) and (2.20) in the HTL approximation, the nonperturbative part of the symmetric propagator in the static limit reads (3.12) The Fourier transform of the propagators leads to the potential 7 8 , Let us make a comment on the derivation. The use of Eq. (3.11) might look problematic, because when we modify the resummed propagator as Eq. (3.9), the self-energies should in general be modified and be different from the perturbative one. The Dyson-Schwinger equation (2.26) can be solved by the following expression, The parameter c can be in fact arbitrary, because the sum of the first two terms on the right hand side is zero. In the thermal equilibrium, the most convenient choice is c = 1 + 2f 0 , where f 0 = f 0 (p 0 ) is the Bose distribution. It is convenient because the FDT holds in the thermal equilibrium, Π S (P ) = (1 + 2f 0 ) sgn(p 0 ) (Π R − Π A ), (3.17) because of which the last two terms of Eq. (3.16) are equal in magnitude. Therefore, in equilibrium, the symmetric propagator can be expressed in three equivalent ways, The one given here is the second model discussed in Ref. [74]. In their first model, the non-perturbative retarded propagator is modeled asD (3.13) The real part of this propagator is m 2 G /(p 2 + m 2 D ) 2 and is the same as the non-perturbative propagator discussed in Ref. [78]. This choice results in the KMS potential [13] for the real part. 8 The same real part of the potential is obtained in Ref. [45] through a different line of reasoning.
If we use the first expression, it is evident that whenD R has an additive contribution as Eq. (3.9), the symmetric propagator gets an additive contribution asD S =D p S +D np S . Namely, Eq. (3.11) does not actually depend on the self-energies.
From this consideration, it might seem difficult to extend this prescription to nonequilibrium. This is because, we do not have the FDT in non-equilibrium, and we have to use Eq. (3.11) to compute the resummed symmetric propagator, but there is no way to determine Π S when there are perturbative and non-perturbative contributions inD R . However, in the current case of the bulk viscous correction, we can determine theD S even though it is non-equilibrium, thanks to the existence of a modified FDT (2.35). We can choose the parameter as c = 2T λ/p 0 and we can expressD S as Because of the modified FDT (2.35), and the last two terms in Eq. (3.19) cancel, just like the case of equilibrium. Therefore, we can express the resummed symmetric propagator in three ways,D (3.20) In the first expression, it is evident that when the retarded (or advanced) propagator gets an additive contribution as Eq. (3.9), the non-perturbative contribution additively contribute toD S =D p S +D np S . In the presence of bulk viscous corrections, let us we consider the same form of the non-perturbative contribution to the retarded self-energy, where Π R is given by Eq. (2.9). Note that the second term does not contribute to the real part of the potential, because its real part vanishes in the static limit p 0 → 0. We can compute the corresponding non-perturbative part of the resummed symmetric propagator using Eq. (3.20) asD np On the left panel of Fig. 5, we plot the imaginary part of the potential for both equilibrium and non-equilibrium cases. Qualitatively, how the imaginary part is affected is the same as the case (3.5). | Im V | gets enhanced at small r, and suppressed at large r.

Approach based on a generalized Gauss law
In Ref. [64], a different prescription to obtain the potential is given. They have used a generalized Gauss law, that can be applicable to a string-like potential, and the HTL permittivity to obtain an analytic form of the in-medium potential. The real part of the potential is the same as Eq. (3.14). The imaginary part of the potential in the thermal equilibrium is given by where G is Meijers G-function. We can use the same prescription straightforwardly to obtain the potential in the presence of the bulk viscous corrections, by using the modified dielectric permittivity (2.31). The real part is again just a replacement m D → m D,R . The imaginary part is given by In this prescription too, the way the bulk correction affects the potential is also very similar to the case (3.5). At small r, the second term is r 4 while the first term is r 2 ln r, so the first term is dominant, where its magnitude is enhanced for Φ > 0. At large r, the second term is dominant. On the right panel of Fig. 5, we plot the imaginary part (3.25) with and without bulk corrections. For Φ > 0, | Im V LR | is enhanced in the small r region, and suppressed at large r.

Comparison of the real part of the potentials
Let us conclude this section by giving a comparison of several models of heavy quark potentials in the thermal equilibrium. In Fig. 6 Figure 6. Comparison of the models of the real part of the in-medium potential.
V TKP ) and (3.14) (denoted by V GDPM/LR ) in the absence of bulk viscous corrections. We also plotted one of the first models of the finite temperature potential is given by [13] V All the models are parametrized by the Debye mass m D , strong coupling constant α, and the string tension σ. We take common values for the plot. Namely, all the models asymptote to the Cornell potential with the same parameters in the small distance limit. For a comparison, we also plotted the perturbative HTL potential (the first term of Eq. (3.4)), the Cornell potential, and Coulomb potential. The potentials Re V KMS , Re V TKP , and Re V GDPM/LR start to deviate from the Cornell potential around the inverse Debye mass and end up between the Cornell and Coulomb potential.

Heavy quarkonia in the presence of bulk viscous corrections
In this section, we discuss the effect of bulk corrections on the properties of heavy quarkonia, such as the binding energies and decay widths, based on the potential obtained in the previous section.

Computational setup
Let us describe the computational procedure. In order to study the in-medium properties of quarkonia, we here solve the Schrödinger equation for a heavy quarkonium to obtain the wave function, using the real part of the in-medium potential. The potential is the function of only the radial coordinate and we only have to solve the ordinary differential equation of the radial part of the wave function. The time-independent Schrödinger equation for the radial wave function reads where m q is the reduced mass of the quarkonium system. We numerically solve this using the real part of the potential (3.4) modified by the bulk viscous correction to obtain the wave functions and eigenvalues 9 . The binding energy is given by the difference between the asymptotic value of the potential and the eigenvalue , In the case of the potential (3.4), the asymptotic value is Using the imaginary part of the potential (3.5), we can make an estimate for the thermal decay width Γ by In a thermal environment, the wave function cannot exist as a steady state, but is transient. This way of the estimate of the decay width treat the imaginary part as a perturbation and its validity worsens when the decay width become comparable the binding energy. Below is how the parameters are set: • For the Debye masses, we use the perturbative expressions (2.10, 2.22) for different temperatures/chemical potentials, including the bulk viscous corrections. Since potentials are parametrized by the Debye mass one could regard the Debye masses as a parameter and fit it to reproduce the potential computed from the lattice QCD data. We do not take this approach here and use the HTL expression, since the primary focus in this study is to understand the nature of bulk viscous modification on heavy-quark properties.
• We set the reduced heavy quark mass to m q = 1.25/2 GeV for cc and m q = 4.66/2 GeV for bb.
• For the coupling constant, we use the one-loop result, 5) with N c = N f = 3. Namely the renormalization scale is taken to be 2π The scale Λ is chosen to Λ = 0.176 GeV requiring α S (1.5 GeV) = 0.326 is satisfied to match the lattice measurements [79]. 9 The interpretation of the role of a complex potential needs some care. As shown in [31], the complex potential dictates the time evolution of unequal time point-split meson-meson correlator, which is related to spectral functions, and the binding energies and decay widths can be read off from the spectral functions. Strictly speaking, the potential does not describe the time evolution of the wave function itself. In this study, we used the real part to dictate the wave function itself and computed the binding energies and decay widths based on the wave function. An advantage of this approach is that it gives us an intuitive picture on the behavior of heavy quarkonia in the medium. The resultant melting temperature is in agreement with those obtained from the approach [31]. When the decay width becomes comparable to the binding energy, it would be better to directly compute the spectral function.

Wave function
Let us first discuss the qualitative feature of the wave functions of quarkonia from the finite temperature potential (3.4). At r → 0, the potential approaches a Coulomb potential with coefficient α, On the other hand, at large distances, the potential also approaches a Coulomb potential, but with a different coefficient, with α ≡ 2σ/ m 2 D,R . The switching of those two regimes happens around r given by the inverse Debye mass. In the limit of the large Debye mass (or high temperature), this part is flattened since the coefficient α goes to zero, and the potential is dominated by the screened HTL one. In the left panel of Fig. 7, we plot the potential as well as its asymptotic Coulomb potentials. In the right panel of Fig. 7, we show the ground state wave function from the potential (3.4) normalized at r = 0 , as well as those from the asymptotic Coulomb potentials, at T = 0.3 GeV. At this temperature, the Debye mass ∼ 0.63 GeV is comparable to the reduced mass of the charm quark, and the wave function is away from both of the two Coulomb wave functions.
In Fig. 8, we show how the wave function is deformed as the temperature goes higher, for the ground state and the first excited states with = 0. A larger temperature results in a larger Debye mass, because of which the wave function is more delocalized and the size of the quarkonium becomes larger. Since the bulk viscous correction comes through the modification of the Debye mass, for the real part of the potential. When Φ > 0, the Debye mass becomes larger, and the wave function becomes more delocalized.  Figure 9 shows the binding energies and decay widths for J/ψ, Υ, and Υ states. As a check of the systematic dependence, we have changed the value of the scale Λ by factors of 1/2 to 2, and different lines The difference in behavior is related to the size of the wave function. In Fig. 10, we show the mean radius of the quarkonia states,r ≡ r 2 , with r 2 = dr r 2 |ψ(r)| 2 r 2 dr r 2 |ψ(r)| 2 . (4.8) At around T = 0.25 GeV, the mean radius of Υ becomes larger than the inverse Debye mass. When the wave function is small compared to the inverse Debye mass, the wave function is not sensitive to the screening and its behavior is determined by the Cornell-like potential. In the limit of tight wave function (large mass), the shape of the wave function is qualitatively close to the Coulomb wave function. Then, the binding energy is given by Since α is an increasing of Λ, the binding energy increases for larger Λ, ∂α ∂Λ > 0. On the other hand, at higher temperature, the size of the wave function grows and the Debye mass becomes relevant. The potential asymptotically approaches a Coulomb potential with a different coefficient α ≡ 2σ/ m 2 D,R . In this situation, the binding energy reads Therefore, when the size of the wave function is comparable or larger than the Debye screening length, the binding energy is a decreasing function of Λ. This is why the Λdependence of the binding energy changes around T = 0.25 GeV for Υ. In the case of J/ψ and Υ , the mean radius is larger than the inverse Debye mass, and the binding energy is always an increasing function of the scale Λ.

Effect of bulk viscous corrections
Now let us discuss the effect of the bulk viscous corrections on the binding energies and decay widths. In Fig. 11, we present the binding energies and decay widths for charmonium (top) and bottomonium (bottom) states as a function of temperature. The solid lines correspond to the cases without bulk correction, and the dotted lines are with the bulk correction. Generically, the bulk correction for Φ > 0 lowers the binding energies, since the bulk viscous correction make the Debye mass heavier and the screening becomes stronger. For a negative Φ, the effect is opposite.
As for the decay widths, as shown in the right column of Fig. 11, the magnitude of the decay width is enhanced in the presence of bulk viscous corrections. The decay width of the ψ state with the bulk correction (top-right) shows flattening at higher temperatures  around 0.4 GeV. This is related to the nature of the imaginary part discussed in Sec. 3.1.
As we discussed there, in the presence of the bulk viscous correction, | Im V | is enhanced at small r, while it is suppressed at higher r. The wave function of J/ψ is smaller in size, and in this region | Im V | is increased, resulting in the larger decay width. On the other Figure 12.
Shift of the melting temperature. The temperature dependence of the binding energies (solid lines) and decay widths (dashed lines) of J/ψ for different parameters Φ (indicated by different colors) is shown. The position of the intersection of a solid and dashed lines of the same color indicate the melting temperature for the corresponding parameter Φ. hand, the excited states are larger and are more sensitive to the large r part of | Im V |.
The enhancement of the decay width of ψ become saturated, because the wave function is now in the region where | Im V | is suppressed by the bulk correction. Based on the computations of the binding energies and decay widths, we can make an estimate of the melting temperature T melt of quarkonium states. We adopt a common criterion that the binding energy coincide with the decay width, E bin (T melt ) = Γ(T melt ). For example, in Fig. 12, we plot the E bin (solid line) and Γ (dashed line) of J/ψ states. Different color corresponds to different parameter Φ. The position at which a solid and dashed lines of the same color intersect is the melting temperature for the parameter Φ. In Table 1, we compare the melting temperatures computed in the current setup (without bulk correction) as well as the results from Ref. [64], that are based on the extraction of the potential with lattice QCD data and subsequent computation of spectral functions. In order to gain a sense of uncertainty, we varied the scale Λ around 0.176 GeV, and the  Table 1.
Melting temperatures for different states in the absence of bulk viscous correction (Φ = 0) , compared with a recent work based on the lattice QCD [64].
upper and lower numbers for our data in the table correspond to Λ = 1/2 × 0.176 GeV and Λ = 2 × 0.176 GeV, respectively. Although there is essentially only one parameter σ (the coupling α S is given by the one-loop result and the Debye mass is given by the HTL) for the case of no bulk viscous correction, the computed melting temperatures is in a reasonable agreement.
In Fig. 13, we plot the melting temperatures of the states J/ψ, Υ and Υ as a function of the bulk viscous parameter Φ. The colored regions indicate the range between the values of Λ, 2×0.176 GeV and 0.5×0.176 GeV, and the solid line is for Λ = 0.176 GeV. We observe a mild decrease of the melting temperature as a function of Φ. The slope of Υ is steeper compared to J/ψ and Υ. Therefore, we have found that the bulk viscous effect indeed affect the melting temperature, especially for excited states. When the plasma is close to the critical point, the bulk viscous contribution is expected to be enlarged. This can lead an anomalous behavior of observables related to heavy quarkonia in the beam energy scan program, although at this state it is difficult to make a quantitative estimate. In order to make a quantitative predictions, one should combine such behavior of quarkonia with a dynamical framework based on hydrodynamics.

Summary and discussions
Heavy quarkonia are useful in probing the nature of the medium around them, through the modification of their properties. In this article, we have studied how the non-equilibrium bulk viscous corrections are imprinted in the properties of heavy quarkonia. The bulk viscous correction modifies the distribution functions of thermal particles, from which the modified dielectric permittivity is computed in the HTL approximation. Using the dielectric permittivity, we computed the heavy quarkonia potential and examined how the potential is deformed in the presence of bulk viscous corrections. As for the real part of the potential, the bulk viscous correction parametrized by Φ > 0 leads to a larger Debye mass, hence stronger screening. The magnitude of the imaginary part is enhanced at short distances, and is suppressed at large r, as shown in Fig. 4. We have tried three prescriptions to accommodate the non-perturbative string-like potential at finite temperatures. Qualitative features of how the bulk viscosity affect the potential is found to be the same among different prescriptions, which indicates the robustness of the results.
We solved the Schrödinger equation with the real part of the potential to obtain the deformed wave functions. We computed the binding energies and decay widths for J/ψ, ψ , Υ, Υ states at different temperatures and the parameter Φ that quantify the bulk viscous corrections. Basically, a positive Φ leads to a Debye mass and it reduces the binding energy. On the other hand, the decay width is enhanced for a given temperature. Because of those, the melting temperature, at which the binding energy equals the decay width, is reduced for Φ > 0 and is enhanced for Φ < 0. When the system is near the critical point, the bulk viscous effect is expected to be enhanced, that would affect the melting temperature. It would be interesting to find an anomalous behavior of observables such as R AA in the beam energy scan program.
Finally, let us make several comments about the possible future directions. A unique ability of potential models is that it can be extended to non-equilibrium situations. It would be interesting to combine the method of potentials with the technique of the QCD sum rule [80][81][82][83][84] to gain insight into the properties of out-of-equilibrium QCD.
Technically and conceptually, in order to understand the dynamical evolution of heavy quarkonia as a quantum state, the analysis based on the theory of open quantum systems would be desirable, which is actively studied recently [85][86][87][88][89][90][91][92]. It would be interesting to apply/derive such a framework to non-equilibrium environments which can have anisotropic noises and see how the time evolution of quarkonia is affected.
Although we have demonstrated that heavy quarkonia have a potential to be sensitive to the bulk viscous nature of the surrounding media, how to experimentally probe this is a nontrivial question. Due to the non-equilibrium correction, the fluctuation-dissipation theorem is violated, which leads to two different Debye masses for retarded and symmetric propagators. If we can make an estimate of those two mass scales independently, it will be possible to explore the non-equilibrium bulk viscous corrections, that are expected to be significant near the critical point.