Spectral functions of heavy quarkonia in a bulk-viscous quark gluon plasma

We study the properties of quarkonia inside a bulk-viscous quark gluon plasma. The non-equilibrium nature of the medium is encoded in the deformed distribution functions of thermal quarks and gluons, with which we compute the dielectric permittivity within the hard thermal loop approximation at one-loop. The modified dielectric permittivity is used to calculate the in-medium heavy quark potential, and using the potential we compute spectral functions, which reflect the physical properties of heavy quarkonia. We discuss how the bulk viscous effect influences quantities such as binding energies and thermal widths. Based on those properties, we discuss the implications of the bulk viscous effect on the physical observables such as ψ′ to J/ψ ratio and the nuclear modification factor, RAA. In particular, we argue that the nuclear modification factors of excited and ground states show different sensitivities to the bulk viscous nature of a plasma, which is potentially useful for the critical point search.


Introduction
The suppression of heavy quarkonia [1,2] is one of the first proposed signals for the formation of a quark gluon plasma (QGP) in heavy ion collision experiments. The properties of heavy quarkonia in the vacuum are well-described by non-relativistic potential models [3] with a Coulombic behaviour at short distances and a non-perturbative linear rise at large distances [4][5][6][7][8][9]. At finite temperatures, inter-quark forces are modified because of the medium effects. It has been argued that above the crossover temperature, T c , the screening is strong enough to dissociate J/ψ state [10]. The melting of quarkonia states is predicted to occur sequentially [11]. The heavy ion collision experiments at RHIC and LHC have also observed relative suppressions of higher states for bottomonia [12,13] and charmonia [14].
In this work, we study the effect of bulk viscosity on the properties of heavy quarkonia and its implications on experimental observables such as the ψ /J/ψ ratio and the nuclear JHEP02(2022)207 modification factor R AA . We incorporate the bulk viscous correction through the modification of quasiparticle distribution functions in the momentum space. This in turn modifies the dielectric permittivity, which can be used [56] to obtain an in-medium potential for a bulk-viscous medium. 1 Using the modified heavy quark complex potential, we compute the spectral functions [64] of quarkonium states in a bulk viscous medium. From those spectral functions, we can read off physical properties such as binding energies, resonance masses and decay widths of quarkonium states. We also compute ψ /J/ψ ratio and nuclear modification factor R AA , which are measurable in heavy ion collision experiments.
Let us state the difference of the current work with our former one [56] in which the effect of bulk viscosity on heavy quarkonia is studied. In ref. [56], the dielectric permittivity is computed in the Hard Thermal Loop (HTL) approximation using the thermal distribution of massless quarks and gluons. Here, we employ a quasiparticle model, in which quasiparticles have masses that depend on the temperature and chemical potential. This model has a better description of thermodynamic quantities near the crossover temperature. In this work, we compute the spectral functions with a complex potential in medium with a method developed in ref. [64], while in ref. [56] we solved the Schrödinger equation using the real part of the potential and used the solved wave function to compute decay widths, which may not be reliable when the imaginary part of the potential is large. Here, we will also discuss implications on experimental observables, ψ /J/ψ ratio and R AA .
The rest of this article is structured as follows. In section 2, we give the derivation of Debye masses and color dielectric permittivity for a bulk viscous medium. In section 3, we discuss the spectral functions of quarkonia in a bulk viscous medium. In section 4, we discuss how the physical observables, ψ /J/ψ ratio and nuclear modification factor R AA , are affected by the bulk viscous corrections. Section 5 is devoted to a summary.

Color dielectric permittivity of a bulk viscous QCD plasma
In this section, we give the computation of the color dielectric permittivity of a bulk viscous medium, which will be used later to get an in-medium potential. First, we introduce the quasiparticle model, and then we discuss how we compute the Debye masses and dielectric permittivity from the modified gluon self energies and propagators following ref. [65].

Quasiparticle model
The quasiparticle model has been used in the study of the thermal properties of a QGP. It reproduces the behaviours of thermodynamic quantities near the cross over temperature T c , where the perturbative QCD is not reliable. In this model, the system of interacting massless partons (quarks and gluons) is effectively described as an ideal gas of massive non-interacting quasiparticles. The quasiparticle mass, m, is given a dependence on the temperature and chemical potential, which arises due to the interactions of partons with the surrounding medium. Such a functional dependence of quasiparticle mass turns out 1 The modified complex potential contains both the perturbative Coulombic as well as the non-perturbative string-like terms. There are different ways for incorporating non-perturbative terms in the potential [20,21,40,[57][58][59][60]. Non-perturbative contributions are also observed in lattice QCD [61][62][63].

Debye masses and color dielectric permittivity
In this subsection, we compute the Debye masses and dielectric permittivity of a bulk viscous plasma. We introduce the bulk viscous correction by deforming the distribution functions of thermal quasiparticles. We parametrize the deformed distribution function in the following way [65] where f id is the isotropic reference distribution (without non-equilibrium corrections) and k ≡ 1 Here +(−) sign is for Bose (Fermi) distribution, m is the quasiparticle mass, whose temperature-and chemical potential-dependence is given by eq. (2.1). Here, the parameter Φ quantifies the strength of the bulk-viscous correction. We here compute the propagators in the presence of the correction (2.3). It will turn out that the effect of the bulk viscous correction is to shift the Debye masses for the retarded, advanced, and symmetric self energies.
Let us first look at the equilibrium contribution to the retarded self energy. In the HTL approximation, the one-loop quarks contribution to the retarded self energy, Π R (P ) in the presence of mass for Φ = 0 is written as [81]

4)
2 This form is motivated by the following deformation of the distribution function [65,80]

JHEP02(2022)207
where f id ± (k) are the distribution functions for quarks/antiquarks in the thermal equilibrium, which is given by The integral can be performed to give (2.7) where Li n (x) denotes the polylogarithm function. Similarly, one can compute the gluon contribution to the retarded self energy using the gluon distribution function in the presence of quasiparticle mass (defined in eq. (2.1)) as where f g (m) is given by (2.10) Therefore, the total retarded self energy at the thermal equilibrium is [65] Π id R (P ) = Π id,q R (P ) + Π id,g R (P ) where m D,R is the retarded Debye mass at Φ = 0, In the massless case, above equation reduces to the familiar expression,

JHEP02(2022)207
Let us now look at the bulk viscous correction. The quark contribution to the retarded self energy in the presence of bulk correction reads (2.14) and the gluon contribution is The bulk viscous correction to the retarded self energy can be written as where δm 2 D,R is given by Hence, the total retarded self energy in the presence of bulk viscous corrections can be written as where m 2 D,R = m 2 D,R + δm 2 D,R . The retarded propagator can be obtained from the retarded self energy. Using eq. (2.18), we can calculate the temporal component of the resummed retarded propagator in the static limit, p 0 → 0, as (2.19) The advanced propagator can be obtained by the complex conjugate of the retarded propagator.
In a similar manner, we can compute the symmetric self energy Π S (P ) and the symmetric propagator D S (P ). The quark contribution to the symmetric self energy for Φ = 0 [65,81] is

JHEP02(2022)207
where Θ(x) is the step function, and K n (x) denotes the modified Bessel function of the second kind. The gluon contribution to the symmetric self energy gives The total symmetric self energy at Φ = 0 is where m D,S is the symmetric Debye mass, In the absence of a quasiparticle mass (m = 0) and a bulk viscous correction (Φ = 0), the symmetric Debye mass reproduces the familiar expression, To the first order in Φ, the bulk viscous correction to the quark contribution of symmetric self energy reads 25) and the correction to the gluon contribution is The total bulk viscous correction to the symmetric self energy is where δm 2 D,S is defined by Thus, the effect of the bulk viscous correction enters as a shift of the Debye mass, and the total symmetric self energy can be written as In the computation of the Debye masses, the coupling constant g 2 = 4πα s is evaluated with the one-loop expression, and we take N c = 3, N f = 3, Λ = 176 MeV and M ≈ 3.7T [67]. We can calculate the resummed symmetric propagator by using the relationD S = D R Π SDA . In the p 0 → 0 limit, we havē (2.31) Figure 1 shows the Debye screening masses, m D,R and m D,S , divided by T as a function of the temperature for different values of Φ. The Debye masses deviate from the linear dependence on T at lower temperatures, where non-perturbative effects are important. The Debye masses are increasing functions of the bulk viscous parameter Φ. In figure 2, we compare the Debye masses computed here with the Debye mass computed via the fitting of the potential with lattice QCD data [40]. The Debye masses m D,R and m D,S show qualitatively similar behaviour as those in ref. [40]. We find that the symmetric Debye mass is slightly larger than retarded Debye mass in the presence of quasiparticle mass. The effect of chemical potential is quantitatively small when µ < T , and hereafter we take µ = 0.
The dielectric permittivity, ε(p), can be computed as [35,56] ε −1 (p) = lim Using eqs. (2.19) and (2.31) into eq. (2.32), the dielectric permittivity is written as (2.34) In the limit of vanishing mass and bulk viscous correction, both the Debye masses m D,R and m D,S approaches to the equilibrium Debye mass m D and the equilibrium expression is reproduced. The dielectric permittivity (2.34) will be used in the computation of the in-medium heavy quark potential.

Heavy quarkonia in a bulk viscous medium
In this section, we discuss the effects of bulk viscous correction on the physical properties of heavy quarkonia. First, we compute the in-medium heavy quark potential with the modified permittivity computed in the previous section. Based on the potential, we compute the spectral functions, from which the physical properties of quarkonia can be extracted.

In-medium heavy quark potential
In this subsection, we discuss how the heavy quark potential is modified in a bulk viscous medium. Following previous approaches [20,56] based on the linear response theory, we compute the in-medium heavy-quark potential by modifying the Cornell potential using the dielectric permittivity ε(p) as potential in the momentum space. The Fourier transform of V (p) is given by The Cornell potential is parametrized as c is the strong coupling constant and σ is the string tension. Here we take, α = 0.513 and σ = (0.412 GeV) 2 [40]. The Fourier transform of eq. (3.2) is given by After substituting eq. (2.34) into eq. (3.1), we obtain an in-medium complex heavy quark potential. Its real part reads In addition to eq. (3.4), we add a constant c, whose value we take c = −0.161 GeV [40]. The first term in eq. (3.4) is the perturbative Coulombic part with a screening, and the second term represents non-perturbative contributions. In the small distance limit (r → 0), eq. (3.4) approaches to the Cornell potential. The bulk viscous correction enters through the modification of the Debye mass, m D,R → m D,R . The imaginary part of the potential is written as where the dimensionless parameter λ is defined as The functions φ n (x) and χ(x) are defined as The function χ(x) is monotonically increasing with χ(0) = 0, and is logarithmically divergent at large x. We regularize this by modifying the following part in eq. (3.8) as with ∆ 3.0369 [40].
In figure 3, we plot the real and imaginary part of the potential. We find that the real part of the potential becomes more flattened for a larger Φ, due to an increased screening. As for the imaginary part, its magnitude is suppressed at larger values of r for Φ > 0.
Let us note the difference in the behaviour of the potential of the current and previous [56] work. The main difference is that we have introduced quasiparticle masses in the current work, while the thermal particles were massless in the previous one. This results in the different temperature dependence of Debye masses. While the Debye masses are always linear in temperature in the previous work, they show non-linear dependence at lower temperatures in the current work (see figure 1) and is suppressed in this region. Thus, the screening effect in the potential is smaller in the present work, which leads to larger values of binding energies as we see later.

In-medium spectral functions of heavy quarkonia
We here discuss the spectral functions of heavy quarkonia in a bulk viscous medium. From the spectral functions, we can extract physical properties of heavy quarkonia such as binding energies and decay widths. For the computation of the spectral functions, we employ the Fourier space method developed in ref. [64]. In this method, the quantity of interest is the unequal-time point-split meson-meson correlator, G > (t; r, r ). The vector-channel spectral function can be obtained from the Fourier transform of the correlator as (ω; r, r ), The correlator is shown to satisfy the following Schrödinger equation, In solving this, we use the potential derived in the previous subsection. Further details of the computation of the spectral functions from eq. (3.12) can be found in ref. [64].
In figure 4, we plot the spectral functions of bottomonium and charmonium for different values of T at Φ = 0. The spectral peaks shift towards lower values of ω at higher temperatures, which is consistent with the results reported in recent works [40,82]. The broadening of peaks occurs because of the imaginary part of the potential, whose magnitude increases at higher temperatures. To quantify the physical properties of each state, we fit each peak of the in-medium spectral functions with the following skewed Breit-Wigner form [40,82,83], where E is the energy of the resonance, Γ is its width, δ is the phase shift, A 1 and A 2 are parameters to account for structures unrelated to the peak of interest. For each state, we obtain the resonance mass and the decay width, and we repeat this process for several different values of T and Φ.  plot the continuum threshold energy (dotted lines), E cont , defined by 2m Q + ReV(r → ∞). The resonance mass of the each state and the continuum threshold energy show a monotonous decrease as a function of T , which is consistent with the previous results [40,82]. The presence of bulk viscous corrections further pushes down the in-medium masses and threshold energy. Figure 10 shows the binding energies of the Υ(1S) and J/ψ(1S) as a function of temperature for different values of Φ. The binding energies are decreasing functions of both T and Φ, which is consistent to our previous results [56]. Quantitatively, the values of the binding energies are higher in the present work, because of the reduced screening effect from the introduction of quasiparticle masses.

���� ���� ���� ���� ���� ���� ���� ���� ���
In contrast, the behaviour of the decay widths is more complicated. In figure 11, we plot the decay widths of the ground and excited states of bottomonium and charmonium for different values of Φ. For the excited states (Υ (2S),Υ (3S), ψ (2S)), the decay widths are decreasing functions of Φ, which is in contrast to our previous computations [56]. However, it does not show a clear tendency for the ground states (Υ(1S), J/ψ(1S)) of both JHEP02(2022)207  bottomonium and charmonium. This behaviour of decay widths can be explained from the following two competing effects: 1 The wave function becomes more spread for larger values of Φ. When |Im V | is small compared to the real part, the decay width Γ is related to the wave function ψ by Γ ∼ − ψ * (ImV )ψ dx, if Φ does not have much effects on ImV , then the spreading of the wave function increases the decay width.
2 Deformation of ImV by Φ. From figure 3, we can see that the magnitude of ImV is suppressed at finite Φ for large r.
The excited states are larger in size compared to the ground states, and are more sensitive to the larger-r part of the |ImV |, and the effect 2 is more important. As a result, the decay widths of the excited states are decreasing functions of Φ. On the other hand, the wave function of the ground states are more compact, and the imaginary part of the potential is not much affected by Φ in this region of r. Hence, none of the two effects is dominant in this case. This behaviour differs from our previous results [56], where the decay widths were increasing functions of Φ for both ground and excited states. The difference arises because of the different behaviour of ImV with Φ at small and large r. In ref. [56], when we have a bulk viscous correction Φ > 0, |ImV | was enhanced at small r and suppressed at large r. The origin of this difference is attributed to the existence of a quasiparticle mass, which was absent in ref. [56].
Thus, we found that, if we look at the decay widths, the excited states are more sensitive to the bulk viscous correction, compared to the ground states. This feature is potentially useful for the search of the critical point, at which the bulk viscous effect is going to be enhanced. It would be interesting to compare the collision energy dependence of R AA of the ground and excited states.
In figure 12, temperatures of the quarkonium states. For a finite Φ, the melting temperature goes down, which is mainly driven by the decrease in the binding energy.

Phenomenological implications
In the previous section, we have discussed how the physical properties of heavy quarkonia are affected by the bulk viscous effect. To discuss their experimental implications, we need to relate them to the observables in heavy ion collisions. In this section, we discuss the effects of bulk viscous correction on physical observables, ψ /J/ψ ratio and R AA .

Relative production yield of ψ to J/ψ
Here, we compute the ratio of the yield of ψ to J/ψ, following the method adopted in refs. [40,60]. The dilepton production rate is written by the spectral function as [84] dR ll

JHEP02(2022)207
where Q q is the electric charge of the heavy quark in units of e, n B is the Bose-Einstein distribution, α e is the fine structure constant, and K = (k 0 , k) is the four momentum. The decaying into dileptons occurs in the vacuum, and we have to carry over the information of the medium to the vacuum spectral function. To achieve this, we perform an instantaneous freeze out at T = T c , under which the in-medium spectral function is projected to the vacuum ones. Let us walk through the procedure. We start with the in-medium lepton-pair production rate (4.1). Since we are interested in the ratio of production yields, prefactors can be dropped, and the yield is proportional to After the change of variables from K to ω(= k 2 0 − k 2 ), the above equation becomes Each peak in ρ V (ω)/ω 2 gives the contribution from each bound state. We fitted each peak structure with the skewed Breit-Wigner (3.14), and determined the resonance mass, M n , and the decay width. We also computed the area A n under each spectral peak. The projection process mentioned above is implemented by replacing the in-medium peaks with delta-function peaks, After performing the integration over ω, the contribution from the state n is written as (4.5) To get the total number density, we also have to divide by the electromagnetic decay rate of a vacuum state, which is proportional to the squared wave function at the origin over the mass squared [96]. The final expression of the ratio is (4.6) We here take M ψ = 3.684 GeV and M J/ψ = 3.0969 GeV, ψ J/ψ (0) = 1.454 GeV 3 , and ψ ψ (0) = 0.927 GeV 3 [97]. Thus, in this procedure, the ratio depends on the resonance mass M n and the area under the peak A n , both of which are affected by medium. It also depends on the temperature at which this "freezing" occurs.
In figure 13, we plot the integrated areas under the bound state peaks for bottomonium and charmonium states. A general trend is that the peak area becomes smaller for larger T and Φ. The proportionality change in the area of the excited state is more pronounced compared to the ground state. Thus, the effect of bulk viscous correction is more significant on the peak areas of the excited states compared to those of the ground states, which is similar to the case of in-medium masses and decay widths.
In order to connect the ratio with the collision energies, we use the freeze-out temperature (which is used in the Bose distribution in eq. (4.5)) fitted to reproduce the particle yields [85], where √ s N N is the numerical value of the center-of-mass energy in GeV. We here set the baryon chemical potential to zero. After substituting all the values in eq. (4.6), we can compute the ψ /J/ψ ratio over a range of center-of-mass energies. Note that our approach applies for the direct production of J/ψ, and the feed down is not considered. Figure 14 shows the relative production yield ψ /J/ψ ratio as a function of center-ofmass energy √ s N N . We compare our results with the statistical hadronization model [85], the Gauss law potential model [40], and experimental data for Pb-Pb [86][87][88][89] and pp [85, 90? -95] collisions. For both of the freezing temperatures, 170 MeV and 200 MeV, our results are close to those from the statistical hadronization model and the Gauss law potential model above 100 GeV center-of-mass energy.
In figure 15, we plot the ψ /J/ψ ratio as a function of Φ. The two lines correspond to different freezing temperatures at which the spectral function is projected to the vacuum one (T = 170 MeV and T = 200 MeV). We find that the ratio moderately increases as a function of Φ at T = 170 MeV, whereas it does not show a clear tendency at T = 200 MeV. The behaviour of ψ /J/ψ ratio as a function Φ can be understood from the behaviour of R ll as a function of Φ for both J/ψ and ψ . From eq. (4.5), we can see that R ll depends on two factors, i.e., areas (A) and masses (M ). The areas and masses decrease as a function of Φ for both the J/ψ and ψ . The decrease in areas as a function of Φ results in the decrease in R ll with Φ. Whereas, the decrease in masses as a function of Φ results in the increase in

JHEP02(2022)207
However, at T = 200 MeV, the effects of Φ on R ψ ll /R J/ψ ll ratio arises due to both the shift in areas as well as the shift in the masses as a function of Φ. Hence, ψ /J/ψ ratio is not showing monotonous behaviour as a function of Φ.

Nuclear modification factor R AA
In this subsection, we give a rough estimate of the nuclear modification factor R AA to see how it gets affected by the bulk viscous correction. The analysis here closely follows ref. [98]. To compute the R AA , we first compute the survival probability of a heavy quarkonium state. We assume that the longitudinal dynamics of the plasma is the Bjorken expansion, and ignore the transverse expansion. We also assume that a produced quarkonium does not move in the transverse plane. We denote the probability that a heavy quarkonium state is in state n by p n . We assume that initially the heavy quarkonium is in state n, p n (τ 0 ) = 1, where τ 0 is the proper time at which the system start interacting with the medium. We treat the interaction of quarkonia with the medium in an adiabatic approximation, on the assumption that the medium evolves slowly. Then, the fate of the state is determined by the decay rate Γ(T (τ )) determined by the local temperature at proper time τ . Under those assumptions, the time evolution of p n (τ ) is given by We choose the initial time as τ 0 = 0.6 fm and the initial condition is p n (τ 0 ) = 1. The survival probability of a given quarkonium state is S = p n (τ f ), where τ f is the final interaction time, which is determined as the time at which the temperature reaches a certain value, T = T f . For an analytical tractability, we fit the decay width (computed in subsection 3.2) with the following function, Γ = aT + bT 2 , (4.9) where a and b are parameters. The survival probability depends on the position in the transverse plane and the impact parameter, S = S (s, b), because of the initial temperature profile. We choose the origin of coordinates in the transverse plane to coincide with the center of one nucleus. The initial energy density, 0 , in the transverse plane is assumed to be proportional to the density of the participants and can be computed by using the optical Glauber model [99][100][101], (4.10) where K is a proportionality coefficient. From the relation between the energy density and the temperature, T ∼ 1/4 , one can obtain the initial temperature profile T 0 (b, s) for a given impact parameter b as

JHEP02(2022)207
where T A is the nuclear thickness function and T 0 (0, 0) is the temperature in the center of the overlapping region for the most central collisions. We consider collisions at the LHC with √ s = 5.02 TeV, and take T 0 (0, 0) = 500 MeV [102,103], A = 207, σ = 70 mb [104]. We use the thickness function of a nucleus with uniform density inside a sphere, where R A = r 0 A 1/3 and r 0 = 1.2 fm. The time dependence of the temperature is given by Using the parametrization of the decay width (4.9), the survival probability can be computed analytically as (4.14) Here, we choose the final temperature to be T f = 158 MeV. For a given impact parameter, we compute the following quantity, The nuclear modification factor R AA for a given impact parameter is computed as , (4.16) where b * is the impact parameter at which the number of participants, N part , is two. We have introduced this normalization so that R AA approaches unity when N part = 2. We compute the nuclear modification factor R AA as a function of N part [105], which can be written as (4.17) Figure 16 shows the computed nuclear modification factor R AA as a function of number of participants N part for Υ and J/ψ. Different lines correspond to different values of Φ. R AA shows a monotonous decrease as a function of N part for both states, as expected. R AA decreases slightly as a function Φ for Υ and increase slightly for J/ψ. We observe that R AA s of the ground states are not much affected by the bulk viscous correction Φ. Since R AA depends on the decay widths and decay widths of the ground states of quarkonia are not much affected by Φ, as shown in figure 11.
However, the decay widths of excited states have a stronger dependence on the bulk viscous correction, as can be seen in figure 11. Thus, we expect that the R AA s of the excited states are going to depend on Φ more strongly than those of the ground states. Within the JHEP02(2022)207  current computational procedure, we cannot make a reliable estimate of them, since we cannot extrapolate the decay widths of the excited states to the temperatures (∼ 500 MeV) realized in the fluid. In order to estimate R AA s of the excited states, a different method is needed. We argue that it will be interesting to compare the collision energy dependence of the R AA s of excited states and ground states to look for the sign of an enhancement of the bulk viscosity near the critical point.

Summary and discussions
In this work, we studied the effect of bulk viscous nature of a plasma on the properties of heavy quarkonia and its implications on the experimental observables in heavy ion collisions. We incorporated the bulk viscous correction, whose strength is quantified by a parameter Φ, by deforming the distribution functions of thermal quarks and gluons. We also introduced quasiparticle masses, with which the behaviour of thermodynamic properties can be well reproduced near the crossover temperature. We computed the dielectric permittivity of a bulk viscous medium within the hard thermal loop approximation at one-loop, and we have obtained an in-medium potential using the dielectric permittivity. From the in-medium potential, we computed the spectral functions, which encode physical properties of quarkonium states.
To quantify the effects of a bulk viscous medium, we have fitted each resonance peak of spectral functions with a skewed Breit-Wigner-type function. By doing this, we can extract quantities such as in-medium masses and decay widths. We found that the in-medium masses and binding energies of the quarkonium states decreases with the increase in bulk viscous correction. The decay widths of the excited states decrease as a function of Φ. On the other hand, the decay widths of ground state do not show a particular tendency as a function of Φ. This can be understood as a result of two competing effects, namely the widening of wave functions because of an increased screening, which tends to increase the decay width, and the deformation of the imaginary part of the potential, which tends to JHEP02(2022)207 decrease the decay width. For excited states, the latter one is dominant, while for ground states the two effects are comparable.
We also studied the implications of the bulk viscous effects on the physical observables such as ψ to J/ψ ratio and the nuclear modification factor R AA . The ψ /J/ψ ratio turned out to show a complicated dependence on the bulk viscous parameter Φ. This appears as a result of the combined effect of the behaviour of the peak area and the resonance mass as a function of Φ. We also gave a rough estimate of the nuclear modification factor R AA for Υ and J/ψ states. We found that R AA s of ground states are not much affected by the bulk viscous correction Φ. However, we expect that the R AA of excited states are more sensitive to Φ, since their decay widths have stronger dependence on Φ than the ground states. Experimentally, it will be interesting to see the collision energy dependence of the R AA s of excited states and ground states. If there is an enhancement of the bulk viscosity near the critical point, that might result in a non-monotonic behaviour of R AA of the excited states.