Effects of event-by-event hydrodynamic fluctuations on bottomonium dynamics in Pb–Pb collisions at $ \sqrt{s_{NN}} = 5.02 $ TeV

High-energy heavy-ion collisions at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) create a deconfined state of strongly interacting matter, known as the quark-gluon plasma (QGP), in which quarks and gluons are no longer confined within hadrons [1]. Heavy quarks and heavy quarkonia are among the most important probes of this medium, because they are predominantly produced in initial hard scatterings and subsequently experience the full space-time evolution of the fireball [2−6]. In particular, bottomonium states provide clean hard probes because the bottom-quark mass is much larger than typical QGP temperatures, making a nonrelativistic description of heavy quarkonium well justified.

From a theoretical perspective, the evolution of heavy quarkonium in hot QCD matter can be formulated in terms of an in-medium heavy-quark potential within effective field theory approaches such as potential nonrelativistic QCD [7, 8]. The corresponding potential becomes complex in the deconfined medium [9−11]. Its real part encodes color screening and modifies the binding structure of quarkonium states, while its imaginary part describes in-medium dissociation induced by interactions with the surrounding light partons. Combined with the time-dependent Schrödinger equation, this framework provides a dynamical and microscopically motivated description of bottomonium suppression in nuclear collisions.

At the same time, relativistic hydrodynamics has achieved remarkable success in describing the soft sector of heavy-ion collisions, indicating that the QGP behaves as an almost perfect fluid with very small specific shear viscosity [12−14]. It is now well established that the medium created in each collision event is not smooth. Instead, event-by-event fluctuations in the initial entropy deposition generate irregular spatial structures and localized hot spots, which are subsequently converted by hydrodynamic expansion into the observed anisotropic flow of final-state hadrons [15−19]. For light hadrons, such fluctuations are essential for understanding higher-order flow harmonics and event-plane correlations.

For heavy quarkonia, however, the quantitative impact of realistic hydrodynamic fluctuations remains less extensively explored. Early studies showed that initial-state fluctuations can noticeably affect the suppression pattern of excited bottomonium states, even when their effect on the ground state remains relatively small [20]. More recently, bottomonium suppression and elliptic flow in fluctuating hydrodynamic backgrounds have been investigated within real-time quantum-evolution frameworks, indicating that the overall influence of fluctuating initial conditions can depend on both the modeling of the medium and the treatment of quarkonium dynamics [21, 22]. These developments make a systematic comparison between bottomonium dynamics in event-by-event fluctuating hydrodynamic backgrounds and smooth hydrodynamic backgrounds timely. In contrast to approaches based on a sharp dissociation temperature or phenomenological inelastic cross sections, the present calculation describes bottomonium suppression through the direct time evolution of the wave function with a complex in-medium heavy-quark potential. In this framework, local temperature variations modify the accumulated real-time evolution and damping of the bottomonium wave packet along its trajectory, rather than removing a state instantaneously once the local temperature exceeds a prescribed threshold. The purpose of the present work is to quantify whether a smooth hydrodynamic background can reproduce bottomonium $ R_{AA} $ and $ v_2 $ obtained from event-by-event fluctuating hydrodynamics within the same time-dependent Schrödinger-equation framework and the same complex-potential uncertainty band. We find that the smooth background reproduces $ R_{\rm{AA}} $ within the theoretical uncertainty, while it systematically underestimates $ v_2 $.

Compared with previous real-time quantum-evolution studies of bottomonium in fluctuating hydrodynamic backgrounds, the present work focuses on a direct comparison between event-by-event fluctuating hydrodynamic media and a smooth optical-Glauber reference background within a Schrödinger-equation framework using a phenomenologically constrained complex heavy-quark potential. We analyze both $ R_{AA} $ and $ v_2 $ for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ across different centrality classes. The remainder of this paper is organized as follows. In Sec. II, we introduce the Schrödinger framework and the in-medium heavy-quark potential. In Sec. III, we describe the event-by-event hydrodynamic background and its validation against soft-hadron data. In Sec. IV, we present the numerical results for bottomonium suppression and flow and discuss their physical implications.

Because of the large heavy-quark mass, the internal dynamics of bottomonium can be described by the time-dependent Schrödinger equation. The radial component of the heavy-quark dipole wave function can be separated as follows:

where $ m_{\mu}=m_b/2 $ is the reduced mass, and $ m_b=4.62\; {\rm{GeV}} $ is the bottom-quark mass. The reduced radial wave function is defined as $ \psi(r,t)=rR(r,t) $, where $ R(r,t) $ is the radial wave function of the quarkonium state. The orbital angular momentum quantum number is denoted by l; $ l=0 $ and $ l=1 $ correspond to the S-wave and P-wave channels, respectively. The time dependence in Eq. (1) enters through the local medium temperature T, which evolves along the trajectory of the propagating $ b\bar b $ dipole in the expanding QGP background.

The in-medium heavy-quark potential is assumed to be complex, $ V(T,r)=V_R(T,r)+V_I(T,r) $. The real part, $ V_R $, describes the screened interaction between the heavy quark and antiquark, whereas the imaginary part, $ V_I $, accounts for in-medium dissociation. For the real part, we use a screened Cornell-type parametrization,

where the Coulomb coupling and string tension are set to $ \alpha=\pi/12 $ and $ \sigma=0.2\,{\rm{GeV}}^2 $, respectively [23]. Here, $ \mu_D $ denotes the Debye screening mass, which is parametrized [24] as $ \mu_D=a_4 T\sqrt{4\pi N_c\left(1+\dfrac{N_f}{6}\right)\dfrac{\alpha}{3}} $, with $ N_c=3 $ and $ N_f=3 $. The parameter $ a_4 $ controls the overall strength of color screening in the medium. For the imaginary part, we adopt the parametrization [24]

where $ \bar r=r/{\rm{fm}} $ is the dimensionless radial distance. $ a_0 $ controls the temperature dependence, whereas $ a_1 $, $ a_2 $, and $ a_3 $ determine the radial dependence of the thermal width. The real part governs the in-medium binding structure, while the imaginary part induces loss of norm of the color-singlet wave packet. In the present work, we explicitly propagate only the color-singlet component of the $ b\bar b $ pair. Color-octet or unbound configurations are not treated as independent degrees of freedom; instead, they are effectively incorporated through the absorptive imaginary part of the potential as probability loss from the singlet sector. The lost components do not contribute to the final projection onto vacuum bottomonium eigenstates. The inverse octet-to-singlet transition, i.e., the regeneration contribution, is not included in the present calculation. The in-medium heavy-quark potential is used in the Schrödinger equation during the QGP phase, defined by $ T>T_c $, with $ T_c=0.17 $ GeV.

In the hadronic gas phase ($ T \lt T_c $), we neglect medium effects and adopt the vacuum Cornell potential in the Schrödinger equation. The parameter ranges of the complex potential are taken from the Bayesian extraction in Ref. [24], namely $ a_0\in[1.001,\,1.407] $, $ a_1\in[0.034,\,0.180] $, $ a_2\in[0.354,\,0.687] $, $ a_3\in[2.000,\,2.579] $, and $ a_4\in[0.100, 0.319] $. These ranges give rise to the uncertainties in $ V_R $ and $ V_I $ shown in Fig. 1, which are used in the subsequent bottomonium calculations. The real component of the heavy-quark potential closely resembles the vacuum Cornell potential, consistent with findings from both Bayesian inference [25] and deep learning approaches [24, 26]. The suppression of quarkonium states is primarily driven by the imaginary part of the potential.

Figure 1.  (color online) The in-medium heavy-quark potential used in the Schrödinger evolution. The left panel shows the screened real part, $ V_R(T,r) $, together with the vacuum Cornell potential, whereas the right panel shows the scaled imaginary part, $ iV_I/T $. The temperature for the bands in $ V_R $ and $ V_I $ is set to $ T=300 $ MeV. The lattice-QCD data are shown for comparison [27].

In relativistic heavy-ion collisions, heavy quarks are produced predominantly through initial hard scatterings. Accordingly, the initial spatial density of $ b\bar b $ dipoles is assumed to be proportional to the density of binary nucleon-nucleon collisions [28].

where $ T_A $ and $ T_B $ denote the thickness functions of the two colliding Pb nuclei, and $ {\boldsymbol{b}} $ is the impact parameter. In the present calculation, this optical-Glauber binary-collision profile is used to sample the initial $ b\bar b $ positions in both the fluctuating and smooth-medium calculations. Thus, the hard-production baseline is kept identical in the two cases, and the comparison focuses on the effect of event-by-event fluctuations in the QGP temperature field. Event-by-event fluctuations of the hard-production profile and their correlations with the MC-Glauber medium fluctuations are not included. The total momentum of each $ b\bar b $ pair is assumed to remain unchanged during the in-medium evolution, so the medium affects only the internal wave function. The transverse momentum of the initially produced $ b\bar b $ dipoles is sampled from a power-law distribution motivated by the measured bottomonium spectra in $ pp $ collisions [25].

where $ \langle p_T^2\rangle=80\; ({\rm{GeV}}/c)^2 $ and $ n=2.5 $.

The time-dependent Schrödinger equation is solved on a radial grid using an implicit finite-difference scheme. At each time step, the discretized radial equation yields a tridiagonal linear system, which is solved using the tridiagonal matrix algorithm [24]. Using the method described above, we evolve each $ b\bar{b} $ dipole on an event-by-event basis along a unique trajectory, initializing the dipoles with random initial positions and total momenta. The Schrödinger evolution continues until the dipole exits the QGP medium at time $ t_f $. The survival fraction of each bottomonium eigenstate $ (n,l) $ is then calculated as $ |c_{nl}(t_f)|^2 $, where the expansion coefficients are obtained by projecting the final reduced radial wave function onto the corresponding vacuum eigenstate, $ c_{nl}(t_f)=\langle \phi_{nl}|\psi(t_f)\rangle $. Here, $ \phi_{nl}(r) $ denotes the radial wave function of the vacuum eigenstate with principal quantum number n and orbital angular momentum quantum number l. The prompt bottomonium yield is obtained by including feed-down contributions from higher excited states [21, 25]. In the present calculation, we evolve the coupled set $ \Upsilon(1S) $, $ \chi_b(1P) $, $ \Upsilon(2S) $, $ \chi_b(2P) $, and $ \Upsilon(3S) $, where the $ \chi_b(nP) $ states denote the spin-averaged triplets. The prompt nuclear modification factor of a final observed state i is calculated as

where $ {\cal{B}}_{j\to i} $ denotes the inclusive feed-down branching fraction from state j to the observed final state i. The feed-down branching fractions are adopted from Ref. [21], and the direct production cross sections listed in Table 1 are obtained by inverting the corresponding feed-down relations.

The nuclear modification factor is calculated using the event-averaged prompt yield as

Here, $ {\mathrm{d}}N^{AA}_{e}/{\mathrm{d}}p_T {\mathrm{d}}\varphi $ denotes the prompt bottomonium yield obtained from a single hydrodynamic event, and $ \langle\cdots\rangle_{\rm{ev}} $ denotes the average over fluctuating hydrodynamic events. For the smooth reference background, the same expression is used without event averaging. The quantity $ {\mathrm{d}}N^{pp}/{\mathrm{d}}p_T {\mathrm{d}}\varphi $ is the corresponding prompt yield in $ pp $ collisions, and $ \langle N_{\rm{coll}}\rangle $ denotes the number of binary nucleon-nucleon collisions in the corresponding centrality class.

The elliptic flow coefficient is calculated with respect to the second-order event-plane angle as

Here, ϕ denotes the azimuthal angle of the bottomonium transverse momentum, and $ \Psi_2^{(e)} $ is the second-order event-plane angle of the e-th hydrodynamic event. After rotating each event to its own event-plane frame, $ \phi'=\phi-\Psi_2^{(e)} $, the angular factor can equivalently be written as $ \cos(2\phi')=\dfrac{p_x^{\prime 2}-p_y^{\prime 2}}{p_x^{\prime 2}+p_y^{\prime 2}} $, where $ p_x' $ and $ p_y' $ are the transverse-momentum components in the event-plane-aligned frame. The elliptic flow coefficient is expected to reflect the anisotropy of the medium distribution arising from event-by-event fluctuations, which is encoded through the interactions of heavy quarkonia with the hot deconfined medium.

Because of the large heavy-quark mass, the internal dynamics of bottomonium can be described by the time-dependent Schrödinger equation. The radial component of the heavy-quark dipole wave function can be separated as follows:

where $ m_{\mu}=m_b/2 $ is the reduced mass, and $ m_b=4.62\; {\rm{GeV}} $ is the bottom-quark mass. The reduced radial wave function is defined as $ \psi(r,t)=rR(r,t) $, where $ R(r,t) $ is the radial wave function of the quarkonium state. The orbital angular momentum quantum number is denoted by l; $ l=0 $ and $ l=1 $ correspond to the S-wave and P-wave channels, respectively. The time dependence in Eq. (1) enters through the local medium temperature T, which evolves along the trajectory of the propagating $ b\bar b $ dipole in the expanding QGP background.

The in-medium heavy-quark potential is assumed to be complex, $ V(T,r)=V_R(T,r)+V_I(T,r) $. The real part, $ V_R $, describes the screened interaction between the heavy quark and antiquark, whereas the imaginary part, $ V_I $, accounts for in-medium dissociation. For the real part, we use a screened Cornell-type parametrization,

where the Coulomb coupling and string tension are set to $ \alpha=\pi/12 $ and $ \sigma=0.2\,{\rm{GeV}}^2 $, respectively [23]. Here, $ \mu_D $ denotes the Debye screening mass, which is parametrized [24] as $ \mu_D=a_4 T\sqrt{4\pi N_c\left(1+\dfrac{N_f}{6}\right)\dfrac{\alpha}{3}} $, with $ N_c=3 $ and $ N_f=3 $. The parameter $ a_4 $ controls the overall strength of color screening in the medium. For the imaginary part, we adopt the parametrization [24]

where $ \bar r=r/{\rm{fm}} $ is the dimensionless radial distance. $ a_0 $ controls the temperature dependence, whereas $ a_1 $, $ a_2 $, and $ a_3 $ determine the radial dependence of the thermal width. The real part governs the in-medium binding structure, while the imaginary part induces loss of norm of the color-singlet wave packet. In the present work, we explicitly propagate only the color-singlet component of the $ b\bar b $ pair. Color-octet or unbound configurations are not treated as independent degrees of freedom; instead, they are effectively incorporated through the absorptive imaginary part of the potential as probability loss from the singlet sector. The lost components do not contribute to the final projection onto vacuum bottomonium eigenstates. The inverse octet-to-singlet transition, i.e., the regeneration contribution, is not included in the present calculation. The in-medium heavy-quark potential is used in the Schrödinger equation during the QGP phase, defined by $ T>T_c $, with $ T_c=0.17 $ GeV.

In the hadronic gas phase ($ T \lt T_c $), we neglect medium effects and adopt the vacuum Cornell potential in the Schrödinger equation. The parameter ranges of the complex potential are taken from the Bayesian extraction in Ref. [24], namely $ a_0\in[1.001,\,1.407] $, $ a_1\in[0.034,\,0.180] $, $ a_2\in[0.354,\,0.687] $, $ a_3\in[2.000,\,2.579] $, and $ a_4\in[0.100, 0.319] $. These ranges give rise to the uncertainties in $ V_R $ and $ V_I $ shown in Fig. 1, which are used in the subsequent bottomonium calculations. The real component of the heavy-quark potential closely resembles the vacuum Cornell potential, consistent with findings from both Bayesian inference [25] and deep learning approaches [24, 26]. The suppression of quarkonium states is primarily driven by the imaginary part of the potential.

Figure 1.  (color online) The in-medium heavy-quark potential used in the Schrödinger evolution. The left panel shows the screened real part, $ V_R(T,r) $, together with the vacuum Cornell potential, whereas the right panel shows the scaled imaginary part, $ iV_I/T $. The temperature for the bands in $ V_R $ and $ V_I $ is set to $ T=300 $ MeV. The lattice-QCD data are shown for comparison [27].

In relativistic heavy-ion collisions, heavy quarks are produced predominantly through initial hard scatterings. Accordingly, the initial spatial density of $ b\bar b $ dipoles is assumed to be proportional to the density of binary nucleon-nucleon collisions [28].

where $ T_A $ and $ T_B $ denote the thickness functions of the two colliding Pb nuclei, and $ {\boldsymbol{b}} $ is the impact parameter. In the present calculation, this optical-Glauber binary-collision profile is used to sample the initial $ b\bar b $ positions in both the fluctuating and smooth-medium calculations. Thus, the hard-production baseline is kept identical in the two cases, and the comparison focuses on the effect of event-by-event fluctuations in the QGP temperature field. Event-by-event fluctuations of the hard-production profile and their correlations with the MC-Glauber medium fluctuations are not included. The total momentum of each $ b\bar b $ pair is assumed to remain unchanged during the in-medium evolution, so the medium affects only the internal wave function. The transverse momentum of the initially produced $ b\bar b $ dipoles is sampled from a power-law distribution motivated by the measured bottomonium spectra in $ pp $ collisions [25].

where $ \langle p_T^2\rangle=80\; ({\rm{GeV}}/c)^2 $ and $ n=2.5 $.

The time-dependent Schrödinger equation is solved on a radial grid using an implicit finite-difference scheme. At each time step, the discretized radial equation yields a tridiagonal linear system, which is solved using the tridiagonal matrix algorithm [24]. Using the method described above, we evolve each $ b\bar{b} $ dipole on an event-by-event basis along a unique trajectory, initializing the dipoles with random initial positions and total momenta. The Schrödinger evolution continues until the dipole exits the QGP medium at time $ t_f $. The survival fraction of each bottomonium eigenstate $ (n,l) $ is then calculated as $ |c_{nl}(t_f)|^2 $, where the expansion coefficients are obtained by projecting the final reduced radial wave function onto the corresponding vacuum eigenstate, $ c_{nl}(t_f)=\langle \phi_{nl}|\psi(t_f)\rangle $. Here, $ \phi_{nl}(r) $ denotes the radial wave function of the vacuum eigenstate with principal quantum number n and orbital angular momentum quantum number l. The prompt bottomonium yield is obtained by including feed-down contributions from higher excited states [21, 25]. In the present calculation, we evolve the coupled set $ \Upsilon(1S) $, $ \chi_b(1P) $, $ \Upsilon(2S) $, $ \chi_b(2P) $, and $ \Upsilon(3S) $, where the $ \chi_b(nP) $ states denote the spin-averaged triplets. The prompt nuclear modification factor of a final observed state i is calculated as

where $ {\cal{B}}_{j\to i} $ denotes the inclusive feed-down branching fraction from state j to the observed final state i. The feed-down branching fractions are adopted from Ref. [21], and the direct production cross sections listed in Table 1 are obtained by inverting the corresponding feed-down relations.

The nuclear modification factor is calculated using the event-averaged prompt yield as

Here, $ {\mathrm{d}}N^{AA}_{e}/{\mathrm{d}}p_T {\mathrm{d}}\varphi $ denotes the prompt bottomonium yield obtained from a single hydrodynamic event, and $ \langle\cdots\rangle_{\rm{ev}} $ denotes the average over fluctuating hydrodynamic events. For the smooth reference background, the same expression is used without event averaging. The quantity $ {\mathrm{d}}N^{pp}/{\mathrm{d}}p_T {\mathrm{d}}\varphi $ is the corresponding prompt yield in $ pp $ collisions, and $ \langle N_{\rm{coll}}\rangle $ denotes the number of binary nucleon-nucleon collisions in the corresponding centrality class.

The elliptic flow coefficient is calculated with respect to the second-order event-plane angle as

Here, ϕ denotes the azimuthal angle of the bottomonium transverse momentum, and $ \Psi_2^{(e)} $ is the second-order event-plane angle of the e-th hydrodynamic event. After rotating each event to its own event-plane frame, $ \phi'=\phi-\Psi_2^{(e)} $, the angular factor can equivalently be written as $ \cos(2\phi')=\dfrac{p_x^{\prime 2}-p_y^{\prime 2}}{p_x^{\prime 2}+p_y^{\prime 2}} $, where $ p_x' $ and $ p_y' $ are the transverse-momentum components in the event-plane-aligned frame. The elliptic flow coefficient is expected to reflect the anisotropy of the medium distribution arising from event-by-event fluctuations, which is encoded through the interactions of heavy quarkonia with the hot deconfined medium.

To simulate the spacetime evolution of the hot QCD medium, we employ the iEBE-VISHNU framework, which provides event-by-event simulations of viscous hydrodynamics for relativistic heavy-ion collisions [32]. In this work, the fluctuating initial conditions are generated using the superMC module based on the Monte Carlo Glauber model. The initial entropy density is constructed as a linear combination of participant and binary-collision contributions,

where $ N_{\rm{part}}(x,y) $ and $ N_{\rm{coll}}(x,y) $ are the local densities of participant nucleons and binary collisions, respectively. The parameter α controls the relative admixture of the two components, and the normalization factor $ s_0 $ is adjusted to reproduce the experimentally measured charged-particle multiplicity.

Event-by-event fluctuations are implemented through stochastic sampling of nucleon positions, together with additional entropy-deposition fluctuations for each participant, modeled by a Gamma distribution in the superMC framework [32]. The Gamma fluctuation parameter is set to 0.75, and the two-component mixing parameter is set to $ \alpha=0.118 $ [33].

The subsequent evolution of the medium is described using a $ (2+1) $-dimensional longitudinally boost-invariant viscous hydrodynamic framework in the Denicol–Niemi–Molnár–Rischke formulation [34], as implemented in the iEBE-VISHNU package [32]. We use a constant specific shear viscosity $ \eta/s=0.08 $, neglect bulk viscosity, initialize the hydrodynamic evolution at $ \tau_0= 0.6\; {\rm{fm}}/c $, and adopt the s95p-PCE equation of state [35]. The hydrodynamic evolution is terminated at the decoupling energy density $ e_{\rm{dec}}=0.18\; {\rm{GeV/fm^3}} $. Particle emission from the decoupling surface is implemented using the Cooper–Frye prescription with shear viscous corrections to the distribution function [36]; the subsequent UrQMD hadronic cascade is not included.

In $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $ Pb–Pb collisions, event-by-event fluctuating backgrounds are generated in several centrality classes and then combined into the centrality intervals used in the bottomonium analysis. In the present study, the bottomonium observables are primarily presented for the 0%−20% and 30%−50% centrality intervals.

Before applying the hydrodynamic background to bottomonium evolution, we briefly verify that the initial setup described in the previous sections provides a reasonable description of representative soft-hadron observables. The final multiplicities of light hadrons predicted by the hydrodynamic model and measured experimentally are listed in Table 2.

The flow coefficients $ v_2 $, $ v_3 $, and $ v_4 $ for identified light hadrons obtained after Cooper–Frye particlization are shown in Fig. 2. In the present calculation, no subsequent UrQMD hadronic cascade is included. Therefore, Fig. 2 should be regarded as a consistency check of the hydrodynamic background rather than a precision calibration of the soft sector. The charged-particle multiplicities and the overall magnitudes of the representative flow coefficients for soft hadrons are reasonably described, which is sufficient for the purpose of the present study: comparing bottomonium evolution in fluctuating and smooth QGP backgrounds within the same hydrodynamic setup.

Figure 2.  (color online) Elliptic flow coefficient $ v_2 $ and higher-order flow coefficients $ v_3 $ and $ v_4 $ for identified light hadrons, $ (\pi^+ + \pi^-)/2 $, $ (K^+ + K^-)/2 $, and $ (p+\bar{p})/2 $, in 30%−40% Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The solid lines show the iEBE-VISHNU calculations, while the symbols indicate the ALICE data [37].

To disentangle the effect of event-by-event fluctuations from that of the average medium evolution, we compare bottomonium observables obtained in two types of hydrodynamic backgrounds: a fluctuating event-by-event background and a smooth reference background. In the former case, each bottomonium trajectory is evolved in an individual hydrodynamic event with its own local temperature inhomogeneities. In the latter case, the evolution is performed in a smooth background corresponding to the same centrality interval. For the smooth reference calculation, the initial entropy-density profile is constructed from an optical-Glauber distribution within the same centrality class. The normalization of the smooth initial condition is chosen to reproduce the same mean final charged-particle multiplicity as in the event-by-event fluctuating hydrodynamic calculation. The subsequent hydrodynamic evolution is performed with the same transport coefficients, equation of state, starting time, and decoupling condition as in the fluctuating calculation.

In a fluctuating hydrodynamic background, the second-order event-plane angle is determined from the initial entropy density according to the standard participant-plane convention used in event-by-event hydrodynamic calculations [16, 19].

where $ \varphi_s $ denotes the spatial azimuthal angle, and $ \langle\cdots\rangle_s $ represents an entropy-density-weighted average in the transverse plane. For each fluctuating event, $ \Psi_2 $ is calculated and used in the expression for bottomonium $ v_2 $.

Figure 3 shows a representative comparison between a smooth hydrodynamic background and an event-by-event fluctuating background in the 0%−20% centrality class. Compared with the smooth reference profile, the fluctuating event contains both localized hot spots and relatively colder regions. Therefore, the effect of event-by-event fluctuations on bottomonium survival is determined by the competition between enhanced dissociation in hotter regions and reduced dissociation in colder regions along the propagation path.

Figure 3.  (color online) Representative temperature contour maps at $ \tau=0.6 {\rm{fm}}/c $ for a smooth hot medium and a fluctuating hot medium in the 0%−20% centrality class of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $.

To simulate the spacetime evolution of the hot QCD medium, we employ the iEBE-VISHNU framework, which provides event-by-event simulations of viscous hydrodynamics for relativistic heavy-ion collisions [32]. In this work, the fluctuating initial conditions are generated using the superMC module based on the Monte Carlo Glauber model. The initial entropy density is constructed as a linear combination of participant and binary-collision contributions,

where $ N_{\rm{part}}(x,y) $ and $ N_{\rm{coll}}(x,y) $ are the local densities of participant nucleons and binary collisions, respectively. The parameter α controls the relative admixture of the two components, and the normalization factor $ s_0 $ is adjusted to reproduce the experimentally measured charged-particle multiplicity.

Event-by-event fluctuations are implemented through stochastic sampling of nucleon positions, together with additional entropy-deposition fluctuations for each participant, modeled by a Gamma distribution in the superMC framework [32]. The Gamma fluctuation parameter is set to 0.75, and the two-component mixing parameter is set to $ \alpha=0.118 $ [33].

The subsequent evolution of the medium is described using a $ (2+1) $-dimensional longitudinally boost-invariant viscous hydrodynamic framework in the Denicol–Niemi–Molnár–Rischke formulation [34], as implemented in the iEBE-VISHNU package [32]. We use a constant specific shear viscosity $ \eta/s=0.08 $, neglect bulk viscosity, initialize the hydrodynamic evolution at $ \tau_0= 0.6\; {\rm{fm}}/c $, and adopt the s95p-PCE equation of state [35]. The hydrodynamic evolution is terminated at the decoupling energy density $ e_{\rm{dec}}=0.18\; {\rm{GeV/fm^3}} $. Particle emission from the decoupling surface is implemented using the Cooper–Frye prescription with shear viscous corrections to the distribution function [36]; the subsequent UrQMD hadronic cascade is not included.

In $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $ Pb–Pb collisions, event-by-event fluctuating backgrounds are generated in several centrality classes and then combined into the centrality intervals used in the bottomonium analysis. In the present study, the bottomonium observables are primarily presented for the 0%−20% and 30%−50% centrality intervals.

Before applying the hydrodynamic background to bottomonium evolution, we briefly verify that the initial setup described in the previous sections provides a reasonable description of representative soft-hadron observables. The final multiplicities of light hadrons predicted by the hydrodynamic model and measured experimentally are listed in Table 2.

The flow coefficients $ v_2 $, $ v_3 $, and $ v_4 $ for identified light hadrons obtained after Cooper–Frye particlization are shown in Fig. 2. In the present calculation, no subsequent UrQMD hadronic cascade is included. Therefore, Fig. 2 should be regarded as a consistency check of the hydrodynamic background rather than a precision calibration of the soft sector. The charged-particle multiplicities and the overall magnitudes of the representative flow coefficients for soft hadrons are reasonably described, which is sufficient for the purpose of the present study: comparing bottomonium evolution in fluctuating and smooth QGP backgrounds within the same hydrodynamic setup.

Figure 2.  (color online) Elliptic flow coefficient $ v_2 $ and higher-order flow coefficients $ v_3 $ and $ v_4 $ for identified light hadrons, $ (\pi^+ + \pi^-)/2 $, $ (K^+ + K^-)/2 $, and $ (p+\bar{p})/2 $, in 30%−40% Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The solid lines show the iEBE-VISHNU calculations, while the symbols indicate the ALICE data [37].

To disentangle the effect of event-by-event fluctuations from that of the average medium evolution, we compare bottomonium observables obtained in two types of hydrodynamic backgrounds: a fluctuating event-by-event background and a smooth reference background. In the former case, each bottomonium trajectory is evolved in an individual hydrodynamic event with its own local temperature inhomogeneities. In the latter case, the evolution is performed in a smooth background corresponding to the same centrality interval. For the smooth reference calculation, the initial entropy-density profile is constructed from an optical-Glauber distribution within the same centrality class. The normalization of the smooth initial condition is chosen to reproduce the same mean final charged-particle multiplicity as in the event-by-event fluctuating hydrodynamic calculation. The subsequent hydrodynamic evolution is performed with the same transport coefficients, equation of state, starting time, and decoupling condition as in the fluctuating calculation.

In a fluctuating hydrodynamic background, the second-order event-plane angle is determined from the initial entropy density according to the standard participant-plane convention used in event-by-event hydrodynamic calculations [16, 19].

where $ \varphi_s $ denotes the spatial azimuthal angle, and $ \langle\cdots\rangle_s $ represents an entropy-density-weighted average in the transverse plane. For each fluctuating event, $ \Psi_2 $ is calculated and used in the expression for bottomonium $ v_2 $.

Figure 3 shows a representative comparison between a smooth hydrodynamic background and an event-by-event fluctuating background in the 0%−20% centrality class. Compared with the smooth reference profile, the fluctuating event contains both localized hot spots and relatively colder regions. Therefore, the effect of event-by-event fluctuations on bottomonium survival is determined by the competition between enhanced dissociation in hotter regions and reduced dissociation in colder regions along the propagation path.

Figure 3.  (color online) Representative temperature contour maps at $ \tau=0.6 {\rm{fm}}/c $ for a smooth hot medium and a fluctuating hot medium in the 0%−20% centrality class of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $.

In this section, we present the effects of event-by-event hydrodynamic fluctuations on bottomonium suppression and elliptic flow in Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. We compare calculations based on fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 0%−20% and 30%−50% centrality intervals. This comparison allows us to examine how the final observables depend on centrality, the binding strength of the bottomonium state, and the interplay between local temperature inhomogeneities and the global geometric anisotropy of the medium. In all result figures, the shaded bands represent the uncertainty propagated from the allowed parameter ranges of the complex heavy-quark potential.

Figure 4 compares the bottomonium $ R_{AA} $ obtained with fluctuating and smooth hydrodynamic backgrounds in the 0%−20% and 30%−50% centrality intervals. The expected hierarchy, $ R_{AA}(1S) \gt R_{AA}(2S) \gt R_{AA}(3S) $, is preserved in both backgrounds, reflecting the different in-medium stabilities of the three states. As shown in the figure, the discrepancy between the two background models remains small in the 0%−20% centrality interval. For both centrality intervals, the effects of event-by-event medium fluctuations on the angle-integrated bottomonium $ R_{\rm{AA}} $ are marginal compared with the smooth reference calculation. The small discrepancy observed for $ \Upsilon(1S) $ is consistent with its stronger binding and reduced sensitivity to local temperature inhomogeneities. Similarly, the $ \Upsilon(2S) $ and $ \Upsilon(3S) $ states exhibit only minor differences between the smooth and fluctuating medium scenarios. Therefore, as reflected in Fig. 4, event-by-event fluctuations in the hot QCD medium do not lead to a quantitatively significant modification of the angle-integrated observable $ R_{AA} $.

Figure 4.  (color online) Comparison of bottomonium $ R_{AA} $ obtained with fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 0%−20% (upper panels) and 30%−50% (lower panels) centrality intervals of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The shaded bands indicate the uncertainties associated with the in-medium heavy-quark potentials. The experimental data points are from the ALICE and CMS measurements [29, 38].

The different p_T dependence from Ref. [20] can be traced partly to differences in the treatment of the hard-production profile. Moreover, because the present calculation does not impose a sharp dissociation temperature, local hot spots do not remove a bottomonium state instantaneously; instead, they modify the cumulative wave-function evolution through the temperature-dependent complex potential. In the present work, the initial $ b\bar b $ positions are sampled from the same optical-Glauber binary-collision profile for both the fluctuating and smooth backgrounds. Therefore, event-by-event correlations between hard-production hot spots and medium hot spots are not included. This reduces the sensitivity of the ensemble-averaged R_AA to local medium fluctuations across the entire p_T range. For weakly bound excited states at low p_T, the suppression is already strong in the smooth background, leaving less room for an additional visible reduction. At higher p_T, although individual trajectories through varying local temperature structures introduce greater statistical variance, the continuous quantum evolution effectively averages over these spatial inhomogeneities. As a result, the central $ R_{\rm{AA}} $ remains close to the smooth reference within the uncertainty band associated with the potential.

The corresponding comparison of $ v_2 $ in the 30%−50% centrality interval is shown in Fig. 5. The figure shows a clear state dependence: the elliptic flow remains very small for $ \Upsilon(1S) $, becomes more visible for $ \Upsilon(2S) $, and is largest for $ \Upsilon(3S) $, following the hierarchy $ v_2(3S) \gtrsim v_2(2S) \gtrsim v_2(1S) $. In contrast to $ R_{\rm{AA}} $, the elliptic flow shows clear sensitivity to event-by-event fluctuations. The fluctuating background yields systematically larger $ v_2 $ than the smooth background for all three states, and the difference becomes more pronounced for the more weakly bound excited states. For $ \Upsilon(1{\rm{S}}) $, the enhancement is small and comparable to the potential-uncertainty band, whereas for $ \Upsilon(2{\rm{S}}) $ and $ \Upsilon(3{\rm{S}}) $, the fluctuating band lies systematically above the smooth band over the intermediate- and high-$ p_{\rm{T}} $ regions, with the separation for $ \Upsilon(3{\rm{S}}) $ exceeding the uncertainty propagated from the in-medium heavy-quark potential.

Figure 5.  (color online) Comparison of $ v_2 $ between fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 30%−50% centrality interval of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The shaded bands represent the uncertainties associated with the in-medium heavy-quark potentials.

This behavior is consistent with the conventional relation between suppression and elliptic flow: the fluctuating and smooth backgrounds give comparable angle-integrated $ R_{\rm{AA}} $, while the fluctuating medium produces larger $ v_2 $. The enhancement can be traced to the anisotropy driving the bottomonium elliptic flow. In the smooth reference, $ v_2 $ is generated by the average, reaction-plane-like geometry of the optical-Glauber profile. In the fluctuating case, each event is analyzed with respect to its own second-order participant plane $ \Psi_2 $, whose eccentricity is enhanced by initial-state fluctuations; the resulting stronger and more irregular azimuthal temperature anisotropy along the trajectories leads to larger differential suppression between the in-plane and out-of-plane directions. Because the more weakly bound excited states are more sensitive to the surrounding medium, this fluctuation-induced enhancement increases from $ \Upsilon(1{\rm{S}}) $ to $ \Upsilon(3{\rm{S}}) $.

In this section, we present the effects of event-by-event hydrodynamic fluctuations on bottomonium suppression and elliptic flow in Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. We compare calculations based on fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 0%−20% and 30%−50% centrality intervals. This comparison allows us to examine how the final observables depend on centrality, the binding strength of the bottomonium state, and the interplay between local temperature inhomogeneities and the global geometric anisotropy of the medium. In all result figures, the shaded bands represent the uncertainty propagated from the allowed parameter ranges of the complex heavy-quark potential.

Figure 4 compares the bottomonium $ R_{AA} $ obtained with fluctuating and smooth hydrodynamic backgrounds in the 0%−20% and 30%−50% centrality intervals. The expected hierarchy, $ R_{AA}(1S) \gt R_{AA}(2S) \gt R_{AA}(3S) $, is preserved in both backgrounds, reflecting the different in-medium stabilities of the three states. As shown in the figure, the discrepancy between the two background models remains small in the 0%−20% centrality interval. For both centrality intervals, the effects of event-by-event medium fluctuations on the angle-integrated bottomonium $ R_{\rm{AA}} $ are marginal compared with the smooth reference calculation. The small discrepancy observed for $ \Upsilon(1S) $ is consistent with its stronger binding and reduced sensitivity to local temperature inhomogeneities. Similarly, the $ \Upsilon(2S) $ and $ \Upsilon(3S) $ states exhibit only minor differences between the smooth and fluctuating medium scenarios. Therefore, as reflected in Fig. 4, event-by-event fluctuations in the hot QCD medium do not lead to a quantitatively significant modification of the angle-integrated observable $ R_{AA} $.

Figure 4.  (color online) Comparison of bottomonium $ R_{AA} $ obtained with fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 0%−20% (upper panels) and 30%−50% (lower panels) centrality intervals of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The shaded bands indicate the uncertainties associated with the in-medium heavy-quark potentials. The experimental data points are from the ALICE and CMS measurements [29, 38].

The different p_T dependence from Ref. [20] can be traced partly to differences in the treatment of the hard-production profile. Moreover, because the present calculation does not impose a sharp dissociation temperature, local hot spots do not remove a bottomonium state instantaneously; instead, they modify the cumulative wave-function evolution through the temperature-dependent complex potential. In the present work, the initial $ b\bar b $ positions are sampled from the same optical-Glauber binary-collision profile for both the fluctuating and smooth backgrounds. Therefore, event-by-event correlations between hard-production hot spots and medium hot spots are not included. This reduces the sensitivity of the ensemble-averaged R_AA to local medium fluctuations across the entire p_T range. For weakly bound excited states at low p_T, the suppression is already strong in the smooth background, leaving less room for an additional visible reduction. At higher p_T, although individual trajectories through varying local temperature structures introduce greater statistical variance, the continuous quantum evolution effectively averages over these spatial inhomogeneities. As a result, the central $ R_{\rm{AA}} $ remains close to the smooth reference within the uncertainty band associated with the potential.

The corresponding comparison of $ v_2 $ in the 30%−50% centrality interval is shown in Fig. 5. The figure shows a clear state dependence: the elliptic flow remains very small for $ \Upsilon(1S) $, becomes more visible for $ \Upsilon(2S) $, and is largest for $ \Upsilon(3S) $, following the hierarchy $ v_2(3S) \gtrsim v_2(2S) \gtrsim v_2(1S) $. In contrast to $ R_{\rm{AA}} $, the elliptic flow shows clear sensitivity to event-by-event fluctuations. The fluctuating background yields systematically larger $ v_2 $ than the smooth background for all three states, and the difference becomes more pronounced for the more weakly bound excited states. For $ \Upsilon(1{\rm{S}}) $, the enhancement is small and comparable to the potential-uncertainty band, whereas for $ \Upsilon(2{\rm{S}}) $ and $ \Upsilon(3{\rm{S}}) $, the fluctuating band lies systematically above the smooth band over the intermediate- and high-$ p_{\rm{T}} $ regions, with the separation for $ \Upsilon(3{\rm{S}}) $ exceeding the uncertainty propagated from the in-medium heavy-quark potential.

Figure 5.  (color online) Comparison of $ v_2 $ between fluctuating and smooth hydrodynamic backgrounds for $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ in the 30%−50% centrality interval of Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02\; {\rm{TeV}} $. The shaded bands represent the uncertainties associated with the in-medium heavy-quark potentials.

This behavior is consistent with the conventional relation between suppression and elliptic flow: the fluctuating and smooth backgrounds give comparable angle-integrated $ R_{\rm{AA}} $, while the fluctuating medium produces larger $ v_2 $. The enhancement can be traced to the anisotropy driving the bottomonium elliptic flow. In the smooth reference, $ v_2 $ is generated by the average, reaction-plane-like geometry of the optical-Glauber profile. In the fluctuating case, each event is analyzed with respect to its own second-order participant plane $ \Psi_2 $, whose eccentricity is enhanced by initial-state fluctuations; the resulting stronger and more irregular azimuthal temperature anisotropy along the trajectories leads to larger differential suppression between the in-plane and out-of-plane directions. Because the more weakly bound excited states are more sensitive to the surrounding medium, this fluctuation-induced enhancement increases from $ \Upsilon(1{\rm{S}}) $ to $ \Upsilon(3{\rm{S}}) $.

In this work, we investigated the influence of event-by-event hydrodynamic fluctuations on bottomonium suppression and elliptic flow in Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02 $ TeV. The internal quantum evolution of the bottomonium states was described using a time-dependent Schrödinger equation incorporating a temperature-dependent complex heavy-quark potential. The QGP background was simulated with the iEBE-VISHNU event-by-event viscous hydrodynamic framework. By comparing fluctuating and smooth hydrodynamic backgrounds for the $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ states, we established a systematic picture of how hydrodynamic fluctuations modify the observables $ R_{AA} $ and $ v_2 $.

Our results indicate that event-by-event fluctuations have only a marginal effect on the bottomonium $ R_{\rm{AA}} $. This reflects the large masses and substantial binding energies of these states, as well as the fact that the angle-integrated yield is largely insensitive to local temperature inhomogeneities once the smooth background is constrained to the same charged-particle multiplicity. By contrast, the elliptic flow $ v_2 $ is systematically enhanced by fluctuations, with a larger enhancement for the more weakly bound excited states. This behavior arises because the participant-plane anisotropy of the fluctuating medium is larger than the average geometry of the smooth background. Bottomonium $ v_2 $ therefore retains discernible sensitivity to event-by-event medium fluctuations, even though $ R_{\rm{AA}} $ is only marginally affected. Higher-order bottomonium flow harmonics and two-particle correlation observables are expected to be more sensitive to event-by-event medium fluctuations than the angle-integrated $ R_{AA} $ and the geometry-dominated $ v_2 $. A dedicated study of these observables requires higher statistics and will be pursued in future work.

In this work, we investigated the influence of event-by-event hydrodynamic fluctuations on bottomonium suppression and elliptic flow in Pb–Pb collisions at $ \sqrt{s_{NN}}=5.02 $ TeV. The internal quantum evolution of the bottomonium states was described using a time-dependent Schrödinger equation incorporating a temperature-dependent complex heavy-quark potential. The QGP background was simulated with the iEBE-VISHNU event-by-event viscous hydrodynamic framework. By comparing fluctuating and smooth hydrodynamic backgrounds for the $ \Upsilon(1S) $, $ \Upsilon(2S) $, and $ \Upsilon(3S) $ states, we established a systematic picture of how hydrodynamic fluctuations modify the observables $ R_{AA} $ and $ v_2 $.

Our results indicate that event-by-event fluctuations have only a marginal effect on the bottomonium $ R_{\rm{AA}} $. This reflects the large masses and substantial binding energies of these states, as well as the fact that the angle-integrated yield is largely insensitive to local temperature inhomogeneities once the smooth background is constrained to the same charged-particle multiplicity. By contrast, the elliptic flow $ v_2 $ is systematically enhanced by fluctuations, with a larger enhancement for the more weakly bound excited states. This behavior arises because the participant-plane anisotropy of the fluctuating medium is larger than the average geometry of the smooth background. Bottomonium $ v_2 $ therefore retains discernible sensitivity to event-by-event medium fluctuations, even though $ R_{\rm{AA}} $ is only marginally affected. Higher-order bottomonium flow harmonics and two-particle correlation observables are expected to be more sensitive to event-by-event medium fluctuations than the angle-integrated $ R_{AA} $ and the geometry-dominated $ v_2 $. A dedicated study of these observables requires higher statistics and will be pursued in future work.