Measurement of inclusive $J/\psi$ production in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV at STAR

A deconfined state of partonic matter, predicted by Quantum Chromodynamics (QCD) and called the Quark-Gluon Plasma (QGP), is produced in ultrarelativistic heavy-ion collisions at the Super Proton Synchrotron (SPS) [1], the Relativistic Heavy Ion Collider (RHIC) [2, 3, 4, 5] and the Large Hadron Collider (LHC) [6]. Understanding the characteristics of the QGP is one of the main research goals of current high-energy heavy-ion experiments. Among the various probes used to study the QGP properties, quarkonia play a unique role. Quarkonia are bound states of heavy quarks and their anti-quarks, whose masses are significantly larger than the QCD scale. Their production yields in heavy-ion collisions are subject to modifications in the presence of the QGP. The hot medium effects on quarkonia include dissociation due to static color screening of the potential between the heavy quark pair [7], dynamical color screening or collisional dissociation caused by interactions with medium constituents [8, 9, 10], and regeneration from deconfined heavy quarks and anti-quarks [11, 12]. Due to the presence of nuclei in the collisions, there are also unavoidable modifications from Cold Nuclear Matter (CNM) effects. These effects arise from changes in parton distribution functions in nuclei compared to those in free nucleons, energy loss in the colliding nuclei, the Cronin effect, nuclear absorption, and other final state effects such as dissociation by co-movers [13, 14, 15, 16, 17].

The $J/\psi$ meson, the most abundantly produced quarkonium state accessible experimentally, has been studied extensively [18, 19, 20, 21, 22]. The inclusive $J/\psi$ yield includes four main sources: $J/\psi$ mesons produced directly during the initial partonic scatterings and via regeneration, those from decays of excited charmonium states, such as $\chi_{c}$ and $\psi$(2S), as well as those from decays of long-lived $b$-hadrons (called non-prompt $J/\psi$). Precision measurements of $J/\psi$ production yields at mid-rapidity have been achieved in Pb+Pb collisions at center-of-mass energies per nucleon-nucleon pair ($\sqrt{s_{\mathrm{NN}}}$) of 17.3 GeV at the SPS [23, 24], 2.76 and 5.02 TeV at the LHC [25, 26], as well as in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 39, 62.4 and 200 GeV at RHIC [27, 28, 29, 30, 31].

Significant suppression is observed at the SPS and RHIC compared to the scaled $J/\psi$ yields in $p$+$p$ collisions of the same energy [23, 24, 27, 28, 29, 30, 31]. The level of suppression is found to be similar, although the temperature and energy density of the produced medium are significantly different due to an order of magnitude difference in $\sqrt{s_{\mathrm{NN}}}$ [27]. The similarity of the suppression is interpreted as the interplay of the energy-dependent dissociation in the QGP and CNM effects, with the regeneration contribution only becoming significant at RHIC. The importance of the regeneration contribution is corroborated by the observation of decreasing $J/\psi$ suppression from the top RHIC energy to LHC energies [25, 32, 33].

There are several theoretical models that employ different dynamic processes, such as $J/\psi$ mesons considered to be dissociated and regenerated continuously during the medium evolution [34, 35, 36, ZHAO2011114] or to be completely melted above the dissociation temperature and then regenerated at the chemical freeze-out [11, 37]. These models can qualitatively describe the experimental measurements. However, the uncertainties in theoretical calculations remain substantial. This is because the underlying mechanisms responsible for the hot medium and CNM effects are not well understood. The model parameters, such as the dissociation temperatures of quarkonium states, the temperature profile of the medium, the heavy quark production cross section and the dynamic evolution of heavy quarks in the medium, are also poorly constrained, which further contributes to the large uncertainties in theoretical models. All these aspects strongly depend on the collision energy, and thus a fine collision energy scan of quarkonium production in heavy-ion collisions can provide stringent constraints on theoretical calculations, which are essential for inferring the properties of the QGP from quarkonium measurements.

The RHIC beam energy scan program, mainly designed for the exploration of the QCD phase structure and to search for the critical end point [38], enables a collision energy scan of $J/\psi$ production between the SPS energy and the top RHIC energy. The STAR and PHENIX experiments have measured $J/\psi$ suppression in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=$ 39 and 62.4 GeV for mid- and forward rapidities, respectively [30, 39], and no significant energy dependence is observed. However, the experimental uncertainties are large due to the limited statistics available. In 2017, the STAR collaboration recorded a significantly larger sample of Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=$ 54.4 GeV compared to those at 39 and 62.4 GeV, with about a factor of six more events. Measurements of $J/\psi$ production at a new energy with improved precision will help to further constrain both the cold and hot medium effects on the $J/\psi$ production in heavy-ion collisions.

In this paper, measurements of the transverse momentum ($p_{\mathrm{T}}$) and centrality dependence of the nuclear modification factor for inclusive $J/\psi$ at mid-rapidity ($\left|y\right|<1$) through the dielectron decay channel in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV are presented. $\langle p^{2}_{\rm{T}}\rangle$ of inclusive $J/\psi$ meson is extracted as a function of the collision centrality. The collision energy dependence of the inclusive $J/\psi$ production is further studied by comparing the new measurements with the published data and theoretical calculations.

The STAR experiment is a general-purpose detector built to study the formation and characteristics of the QGP. It covers the entire azimuth within the pseudorapidity interval of $\left|\eta\right|<$ 1. The dataset utilized in this analysis contains $1.3\times 10^{9}$ minimum-bias events for Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV collected by the STAR experiment in 2017. The minimum-bias trigger requires coincident signals in the two Zero-Degree Calorimeters (ZDCs) [40, 41], or the two Vertex Position Detectors (VPDs) [42]. Charged particles produced in the collisions are tracked and reconstructed in the Time Projection Chamber (TPC) [43], which is immersed in a solenoidal magnetic field of 0.5 T. The event vertex determined using the tracks in the TPC is required to be within $\pm$60 cm (2 cm) from the nominal center of the STAR detector along the beam (radial) direction. To reject pileup events, the reconstructed vertex position along the beam direction is further required to be within 3 cm of the vertex position reconstructed based on coincident signals in the VPDs, which are fast detectors and thus resilient to pileup effects. To avoid trigger detector inefficiency for events with low multiplicities, only data from the 0-60% most central collisions are analyzed. The total number of events used after event selection is 593 million.

The collision centrality is determined by matching the multiplicity distribution of charged tracks ($dN_{\mathrm{ch}}/d\eta$) in data to that from the Monte Carlo Glauber model [44]. The $dN_{\mathrm{ch}}/d\eta$ is obtained by counting the number of charged tracks within $\left|\eta\right|<$ 0.5 and is corrected for vertex position and luminosity dependences, but not for the tracking efficiency. The average number of binary collisions ($\langle N_{\mathrm{coll}}\rangle$) and number of participants ($\langle N_{\mathrm{part}}\rangle$) for each centrality bin are determined using the Glauber model [45].

The $J/\psi\to e^{+}e^{-}$ decay channel is used to reconstruct $J/\psi$ candidates in this analysis. To assure high quality, accepted tracks in the TPC are required to have at least 20 space points used for their reconstruction. The tracks’ distance of closest approach to the event vertex needs to be less than 1 cm. Furthermore, the ratio of the number of TPC space points used for track reconstruction to the maximum possible number along the track trajectory should be larger than 0.52 to remove split tracks.

The main detectors used for electron identification are the TPC, the Time-of-Flight (TOF) detector [46], and the Barrel Electromagnetic Calorimeter (BEMC) [47]. Unless specified otherwise, the term “electron” denotes both electrons and positrons in the following. The TPC provides particle identification (PID) through measurement of the specific energy loss ($\langle dE/dx\rangle$). The variable used is defined as:

where $(dE/dx)_{\rm{measured}}$ is the measured energy loss, $(dE/dx)^{\mathrm{e}}_{\rm{theory}}$ is the predicted energy loss for an electron based on the Bichsel formalism [48] and $R_{\ln(dE/dx)}$ is the resolution of $\ln(dE/dx)$. Consequently, the $n\sigma_{\mathrm{e}}$ distribution for electrons is expected to follow a Gaussian distribution with a mean of zero and width of one. To ensure good resolution, the number of TPC space points used to calculate $dE/dx$ is required to be no less than 15. In combination with the path length measured in the TPC, the TOF provides the velocity information ($\beta=v/c$) by measuring a particle’s flight time, where $c$ is the speed of light. Electrons can be separated from hadrons due to their larger velocities at a given momentum. Furthermore, the BEMC is used to suppress hadron contamination via the measurement of deposited energies by incident particles and comparison to their momenta measured in the TPC. A TPC track is matched to a BEMC cluster geometrically by projecting to the BEMC surface, and the variable $E_{0}/p$, where $p$ is the particle momentum measured in the TPC and $E_{0}$ is the energy deposition of the most energetic tower in the matched BEMC cluster, is used for PID [49]. For electrons, the $E_{0}/p$ is expected to be peaked around one, while for hadrons it is considerably smaller than one.

Table 1 lists the combinations of detectors and specific cut values used for electron identification in different $p_{\rm{T}}$ ranges. For tracks with $p_{\mathrm{T}}\leq 1$ GeV/$c$, a momentum-dependent $n\sigma_{\mathrm{e}}$ cut is applied, along with the requirement $\left|1/\beta-1\right|<0.025$. For $p_{\mathrm{T}}>1$ GeV/$c$, different combinations of PID selections are used, depending on the availability of TOF and BEMC information, in order to maximize the hadron rejection power while maintaining a high identification efficiency. When the track is only matched to TOF, the selection criteria of $\left|1/\beta-1\right|<0.025$ and $-0.75<n\sigma_{\mathrm{e}}<2$ are applied utilizing both TOF and TPC information. If a track leaves a signal in the BEMC but not in the TOF, the BEMC cut of $0.5<E_{0}/p<1.5$ is utilized in conjunction with the TPC cut of $-1<n\sigma_{\mathrm{e}}<2$. The third scenario is that the PID information from TPC, TOF and BEMC can all be used, and the corresponding electron selection cuts are: $\left|1/\beta-1\right|<0.025$, $-1.5<n\sigma_{\mathrm{e}}<2$ and $0.5<E_{0}/p<1.5$. It is worth noting that the $n\sigma_{\mathrm{e}}$ cut is asymmetric with a tighter and lower boundary on the negative side. This is because charged pions are the most abundant hadron species that contaminate the electron sample and their $n\sigma_{\mathrm{e}}$ values are negative. The lower boundary of the $n\sigma_{\mathrm{e}}$ cut also changes with $p_{\mathrm{T}}$ and the availability of TOF and BEMC information, while the upper boundary is kept at 2.

Figure 1 shows the invariant mass ($M_{\mathrm{ee}}$) distribution of $e^{+}e^{-}$ pairs in 0-60$\%$ Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV. A $p_{\rm{T}}>0.2$ GeV/$c$ cut is applied to exclude coherent photon-induced $J/\psi$ production [50, 51, 52].

The raw yield of the $J/\psi$ meson is obtained by counting the $e^{+}e^{-}$ pairs in the mass range of 2.9-3.2 GeV/$c^{2}$, where a clear $J/\psi$ peak is seen around $M_{\mathrm{ee}}=3.09$ GeV/$c^{2}$, and subtracting the combinatorial and residual background within the mass window. The residual background, originating from correlated heavy-flavor decays and Drell-Yan processes, is estimated by fitting the same-event unlike-sign distribution using the unbinned maximum likelihood method with the following components: i) the $J/\psi$ signal shape described by a Crystal-Ball function [53, 54] whose width is constrained by detector simulations tuned to reproduce realistic momentum resolution; ii) the combinatorial background estimated using the normalized unlike-sign $e^{+}e^{-}$ pairs from mixed events [55]; iii) a first-order polynomial function describing the residual background. The mixed-event unlike-sign distribution is normalized according to the ratio of like-sign distributions between same and mixed events. The fitting result is shown as the red solid line in the upper panel of Fig. 1. Also shown in Fig. 1 is the combinatorial background subtracted unlike-sign distribution, overlaid with the fit results of the $J/\psi$ signal plus the residual background, in the bottom panel. The resulting significance of the $J/\psi$ raw signal, calculated as $S/\sqrt{S+B}$ where $S$ and $B$ are signal and background counts, is 32.8 for $p_{\rm{T}}>0.2$ GeV/$c$ and 0-60% centrality class. The efficiency of counting raw $J/\psi$ yields within the mass range of 2.9 - 3.2 GeV/$c^{2}$ is obtained based on the fitted Crystal-Ball function.

The electron detection efficiency consists of the TPC tracking efficiency, TOF and BEMC matching efficiencies and PID cut efficiencies. The embedding technique is used to evaluate the TPC tracking efficiency as well as the BEMC matching and PID efficiencies [44]. A sub-sample of minimum-bias events, which cover the entire period of data taking, is processed using the standard STAR reconstruction procedure after being embedded with the simulated electron tracks propagated through the GEANT3 simulation [44] of the STAR detector. On the other hand, the TOF matching efficiency as well as the TPC and TOF PID cut efficiencies are evaluated using a pure electron sample in data, originating from photon conversions and light meson Dalitz decays [56]. For an electron track with $p_{\rm{T}}$ around 1.5 GeV/$c$, the TPC track efficiency is about 83%, the TOF matching efficiency is about 7% and the BEMC matching efficiency is about 86% in 0-60% centrality. The $J/\psi$ reconstruction efficiency is determined by convoluting the electron detection efficiency as a function of $p_{\rm{T}}$, $\eta$, and $\phi$ with the $J/\psi$ decay kinematics.

Systematic uncertainties include those in signal extraction, TPC tracking efficiency, TOF and BEMC matching efficiency estimations and electron identification. For the signal extraction, the invariant mass range is varied when obtaining the normalization factor for the unlike-sign distribution in mixed events. The raw $J/\psi$ counts are obtained directly from the fit function instead of the default bin-counting method. Additionally, different fitting ranges for estimating the background contribution are tried. The shape of the residual background is changed from the default first-order polynomial function to zeroth-order or second-order function. Finally, the width of the Crystal-Ball function is varied taking into account the uncertainties in tuning the simulation to match momentum resolution in data. The maximum deviation of all the variations to the default is taken as the uncertainty. For the TPC tracking efficiency, track quality cuts are varied simultaneously in data analysis and in extracting the tracking efficiency from embedding to evaluate the uncertainties. An additional $5\%$ systematic uncertainty is included for each single electron tracking efficiency to account for discrepancies arising from imperfect simulation when compared to the data [57]. In terms of the electron identification, the uncertainties on $n\sigma_{\mathrm{e}}$ and $E_{0}/p$ cut efficiencies are estimated by changing these cuts and checking the impact on the corrected $J/\psi$ yields. The default 1/$\beta$ cut efficiency is calculated via counting the fraction of the electron 1/$\beta$ distribution within $\left|1/\beta-1\right|<0.025$. The difference between the default method and that based on a Gaussian fit to the 1/$\beta$ distribution is taken as the uncertainty. The uncertainty in the TOF matching efficiency is evaluated by comparing efficiency differences when using different cuts to select the electron sample. On the other hand, for the BEMC matching efficiency, the uncertainty originates from the difference between utilizing detector simulation and data-driven approach. All the aforementioned uncertainties are listed in Table 2. Global uncertainties, shown as bands at unity in the figures, are from the uncertainties in $\left\langle N_{\rm{coll}}\right\rangle$ and the $J/\psi$ production cross section in $p$+$p$ collisions.

The efficiency and acceptance corrected invariant yields of inclusive $J/\psi$ within $|y|<1$ in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV are shown in Fig. 2 as a function of $p_{\rm{T}}$ for different centrality classes. The invariant yield is extracted as:

where $B_{J/\psi\to e^{+}e^{-}}$ is the branching ratio of $J/\psi$ decaying into an $e^{+}e^{-}$ pair, $\Delta N_{J/\psi}$ is the raw $J/\psi$ count in each $p_{\rm{T}}$ bin for the considered centrality class, $p_{\rm{T}}$ is the bin center, $\Delta p_{\rm{T}}$ is the $p_{\rm{T}}$ bin width, $\Delta y$ is the rapidity coverage, from -1 to 1, $N_{\rm{MB}}$ stands for the number of MB events in the considered centrality class used in this analysis, and $\epsilon_{total}$ is the total $J/\psi$ efficiency.

Data points are placed at $p_{\rm{T}}$ values whose yields are equal to the average yields of the bins [58]. The $p_{\rm{T}}$ positions are determined by fitting the differential yields iteratively with an empirical function:

where $A$, $B$, and $C$ are free parameters. The horizontal bar around each data point indicates the bin width.

The nuclear modification factor ($R_{\rm{AA}}$) is used to quantify the modification to the $J/\psi$ production and is defined as:

where $d^{2}N_{\rm{AA}}/dp_{\rm{T}}dy$ is the $J/\psi$ yield in A+A collisions and $d^{2}\sigma_{pp}/dp_{\rm{T}}dy$ is the $J/\psi$ cross section in $p$+$p$ collisions. $\sigma^{\rm{inelastic}}_{pp}$ is the inelastic $p$+$p$ cross section. Since there are no experimental measurements for the inclusive $J/\psi$ production cross section in $p$+$p$ collisions at $\sqrt{s}=54.4$ GeV, a data-driven method is used to derive it based on world-wide experimental data from $p$+$p$ to $p$+A collisions ranging between $\sqrt{s_{\mathrm{NN}}}=6.8-7000$ GeV [59]. Different formulas are used to fit the collision energy, rapidity and $p_{\rm{T}}$ dependence of the inclusive $J/\psi$ production cross section, and the difference in the interpolated values at 54.4 GeV between these formulas is taken as the systematic uncertainty for the $p$+$p$ baseline.

The $p_{\rm{T}}$-integrated $R_{\rm{AA}}$ of inclusive $J/\psi$ as a function of $\langle N_{\rm{part}}\rangle$ is shown in Fig. 3 for Au+Au collisions at collision energies of 39, 54.4, 62.4, 200 and 5020 GeV [30, 31, 26]. The newly measured $J/\psi$ $R_{\rm{AA}}$ at 54.4 GeV is consistent with previous results at RHIC within uncertainties, indicating no significant collision energy dependence of $R_{\rm{AA}}$ up to 200 GeV, while the precision is significantly improved. Taking into account statistical and systematic uncertainties from Au+Au measurements, the significance of the deviation in $R_{\rm{AA}}$ between 54.4 and 62.4 GeV varies between 1.6$\sigma$ and 2.5$\sigma$ depending on $\langle N_{\rm{part}}\rangle$. The distributions also hint at an increasing suppression from peripheral to central collisions at RHIC energies, consistent with increasing hot medium effects. On the other hand, the $J/\psi$ suppression at the LHC seems to decrease towards central collisions and is larger than at RHIC energies, which is attributed to a larger regeneration contribution at higher energies and more central collisions. Transport model calculations [60] from the Tsinghua group for 39, 54.4 and 62.4 GeV are shown as dashed curves in Fig. 3, which predict very little difference among different energies as observed in data.

Unlike $R_{\rm{AA}}$, another nuclear modification factor, $R_{\rm{CP}}$, can be used to reflect the relative suppression between central and peripheral collisions. The $R_{\rm{CP}}$ is defined as a ratio of the $J/\psi$ yield in central collisions over that in peripheral collisions:

where $dN_{\rm{AA}}/dy$ is the $p_{\rm{T}}$-integrated $J/\psi$ yield in a certain centrality class. One advantage of using $R_{\rm{CP}}$ is that it does not use the $J/\psi$ yield from $p$+$p$ collisions as the reference, which can only be interpolated from parameterization of world data for 39, 54.4 and 62.4 GeV. Furthermore, correlated systematic uncertainties in $J/\psi$ yield measurement largely cancel in the $R_{\rm{CP}}$. In this analysis, the peripheral 40-60% centrality bin is used as the reference. Figure 4 shows the centrality dependence of $R_{\rm{CP}}$ in Au+Au collisions at RHIC. The boxes on the data points represent systematic uncertainties, mainly from signal extraction, while uncertainties from tracking and electron identification mostly cancel. The remaining uncertainties in $\langle N_{\rm{coll}}\rangle$ after cancellation are indicated as shaded bands at unity. A relative suppression of $J/\psi$ yield in central collisions compared to that in peripheral collisions is observed at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV. No significant collision energy dependence of $R_{\rm{CP}}$ has been observed at RHIC energies within uncertainties.

Figure 5 shows the collision energy dependence of $J/\psi$ $R_{\rm{AA}}$ in central heavy-ion collisions of different species.

No significant energy dependence is seen within uncertainties between 17.3 and 200 GeV in central collisions. Two transport model calculations of the collision energy dependence of inclusive $J/\psi$ $R_{\rm{AA}}$, displayed by the same line style, are shown on the left panel from the Tsinghua group [60] and the right panel from the TAMU group [36]. Blue dash-dot-tripled lines represent the suppression of primordial production due to CNM effects and dissociation in the QGP medium, while the red long dash lines denote the regeneration contribution. The final $J/\psi$ $R_{\rm{AA}}$ calculations from the two groups, considering both the suppression and regeneration, are shown as the black dash lines in the two panels of Fig. 5. The theoretical calculations, starting from 39 GeV, are consistent with the observed energy dependence of $J/\psi$ $R_{\rm{AA}}$, indicating that the $J/\psi$ production in high-energy heavy-ion collisions is an interplay of dissociation in the QGP medium, regeneration, and CNM effects. At 17.3 GeV, the transport model calculation from the Tsinghua group underestimates the experimental measurements.

Figure 6 shows $J/\psi$ $R_{\rm{AA}}$ as a function of $p_{\rm{T}}$ for the 0-60% centrality class at different energies at RHIC [30, 31] and for different centrality classes at 54.4 GeV. The $p_{\rm{T}}$ dependence seems to be flatter at $\sqrt{s_{\mathrm{NN}}}=200$ GeV compared to lower energies, which could be due to a larger regeneration contribution at 200 GeV and larger nuclear absorption at lower energies at low $p_{\rm{T}}$. The $p_{\rm{T}}$ dependence of $R_{\rm{AA}}$ also exhibits a similar distribution at 39, 54.4, and 62.4 GeV within uncertainties. In Fig. 6(b), a slightly larger yield suppression is observed towards central collisions compared to that for peripheral collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV, while the $p_{\rm{T}}$ dependence is similar across different centrality classes. Transport calculations from the Tsinghua group, shown for comparison, are seen to overestimate $R_{\mathrm{AA}}$ below 2 GeV/$c$ for 20-40% and 40-60% centrality classes at 54.4 GeV.

The production of $J/\psi$ in high-energy heavy-ion collisions is an interplay of different hot and cold medium effects, while the $J/\psi$ $p_{\rm{T}}$ spectra shapes are crucial for studying these effects individually since different effects may dominate in different $p_{\rm{T}}$ regions. The second moment of the $p_{\rm{T}}$ distribution ($\langle p^{2}_{\rm{T}}\rangle$) can facilitate the comparison of the measured $J/\psi$ $p_{\rm{T}}$ spectra shapes in different centrality bins at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV as well as at different collision energies. The resulting $\langle p^{2}_{\rm{T}}\rangle$ in different centrality classes of Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV are displayed in Table 3, and no significant centrality dependence is seen.

The $\langle p^{2}_{\rm{T}}\rangle$ in $p$+$p$ collisions at $\sqrt{s}=54.4$ GeV is derived to be $2.53\pm 0.05$ $(\rm{GeV}/c)^{2}$ [59], where 0.05 is the total uncertainty, lower than in Au+Au collisions. This is consistent with the observed $p_{\rm{T}}$ dependence of $R_{\rm{AA}}$, which increases towards higher $p_{\rm{T}}$ as shown in Fig. 6. To quantify the change in $\langle p^{2}_{\rm{T}}\rangle$, $r_{\rm{AA}}$ is used. It is defined as the ratio between $\langle p^{2}_{\rm{T}}\rangle$_AA and $\langle p^{2}_{\rm{T}}\rangle$_pp, and is shown as a function of $\langle N_{\rm{part}}\rangle$ in Fig. 7 for heavy-ion collisions with the collision energy ranging from 17.3 GeV to 5.02 TeV [32, 27, 61, 62, 63, 64]. The $r_{\rm{AA}}$ at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV shows a flat distribution against centrality and follows the trend of the collision energy dependence. On the other hand, $r_{\rm{AA}}$ decreases towards central collisions at 5.02 TeV, due to increased regeneration contribution at low $p_{\rm{T}}$.

We report on the measurements of inclusive $J/\psi$ production at mid-rapidity ($|y|<1$) in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV by the STAR experiment at RHIC. The newly measured $J/\psi$ $R_{\rm{AA}}$ at 54.4 GeV is compatible with previous measurements at 39 and 62.4 GeV with improved precision. A significant suppression of the $J/\psi$ yield is observed in the 0-60% centrality class compared to the $p$+$p$ baseline. Hints of larger suppression towards central collisions and low $p_{\rm{T}}$ are seen, and no significant energy dependence of $J/\psi$ $R_{\rm{AA}}$ in central collisions from 17.3 to 200 GeV is observed within uncertainties. The $\langle N_{\rm{part}}\rangle$ dependence of $r_{\rm{AA}}$ at $\sqrt{s_{\mathrm{NN}}}=54.4$ GeV shows a flat distribution and follows the trend of the collision energy dependence. Transport model calculations incorporating dissociation and regeneration in the QGP and the CNM effects can qualitatively describe the data. These results will provide additional constraints on understanding the interplay between different hot and cold nuclear matter effects, which are essential for probing the QGP properties with quarkonia.

We thank the RHIC Operations Group and SDCC at BNL, the NERSC Center at LBNL, and the Open Science Grid consortium for providing resources and support. This work was supported in part by the Office of Nuclear Physics within the U.S. DOE Office of Science, the U.S. National Science Foundation, National Natural Science Foundation of China, Chinese Academy of Science, the Ministry of Science and Technology of China and the Chinese Ministry of Education, NSTC Taipei, the National Research Foundation of Korea, Czech Science Foundation and Ministry of Education, Youth and Sports of the Czech Republic, Hungarian National Research, Development and Innovation Office, New National Excellency Programme of the Hungarian Ministry of Human Capacities, Department of Atomic Energy and Department of Science and Technology of the Government of India, the National Science Centre and WUT ID-UB of Poland, the Ministry of Science, Education and Sports of the Republic of Croatia, German Bundesministerium für Bildung, Wissenschaft, Forschung and Technologie (BMBF), Helmholtz Association, Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) and Agencia Nacional de Investigación y Desarrollo (ANID) of Chile.