Simulations of lattice quantum chromodynamics(QCD) predict that when the temperature of the system reaches up to about $ 155 $ MeV^[1-3], the quark and gluon degrees of freedom that are usually confined in normal nuclear matter are liberated, and a new state of matter called quark-gluon plasma (QGP) can be formed. Such hot and dense deconfined QGP matter has existed in early Universe a few microseconds after the Big Bang, when the temperature of the Universe is so high that quarks and gluons cannot combine to form hadrons due to their violent thermal motions. QGP can also be created in the laboratories via colliding heavy nuclei at ultra-relativistic energies, such as those performed at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). The primary goal of high-energy heavy-ion collisions is to probe the novel properties of QGP and explore the phase structure of strong-interaction QCD matter at various conditions of temperatures and/or densities. Comprehensive experimental and theoretical studies have shown that the highest temperature achieved at RHIC and the LHC is about 350~480 MeV^[4-8], which is well above the phase transition temperature ($ \sim 155 $ MeV) between normal nuclear matter and deconfined QGP.

Relativistic heavy-ion collision is a dynamically-evolving many-body process. It starts from the collision of two heavy ions which have been accelerated to ultra-relativistic energies. The energetic ion collision deposits a lot of energy and produces a large amount of entropy and particles. After a short pre-equilibrium stage, the system arrives at a nearly thermal equilibrium state, i.e., the hot and dense QGP. The fireball then expands and cools until the temperature reaches down to the phase transition temperature when the system hadronizes into hadrons. The produced hadrons continue to interact with each other until the kinetic freeze-out. Finally, all produced particles fly to the detectors. One important task of relativistic heavy-ion collisions is, by studying the observed distributions and correlations of produced particles, to infer the evolution dynamics of the collisions and probe the macroscopic and microscopic properties of QGP.

In this review, we discuss two most important probes of QGP, soft and hard probes, in relativistic heavy-ion collisions. Soft probes refer to the collective properties of soft bulk matter, which consists of more than 99% of particles produced in collisions. Here we focus on anisotropic collective flows of the QGP. Hard probes mainly refer to the particles carrying large transverse momenta ($ p^{}_{\rm{T}} $), which are usually originated from early hard collisions. Here we focus on high $ p^{}_{\rm{T}} $ and heavy flavor hadrons and jets. We will also discuss the collectivity in small collision systems.

One of the most important evidences of QGP formation in relativistic heavy-ion collisions is the strong collective flow^[9-12], which is anisotropic in the plane transverse to the beam. In a typical nucleus-nucleus collision, the produced QGP matter is anisotropic in the transverse plane due to the intrinsic collision geometry and the fluctuations of the initial states^[13-25]. The initial geometric anisotropy can be quantified by the so-called eccentricity $ \epsilon_n $, defined as follows^[20-22]:

where $ m $=3 for $ n $=1 and $ m\!=\!n $ for $ n \geqslant 2 $. $ \Psi_n $ is the orientation of the $ n $-order eccentricity $ \epsilon_n $, $ \{ \cdots \} $ represents the average over the entropy (energy) density profile, and $ (r_\perp, \varphi) $ represents the polar coordinates. Due to the strong interaction among QGP constituents, the initial geometric anisotropy can be translated into the final state anisotropy of particle momentum distribution. The momentum anisotropy can be quantified via the Fourier decomposition of the spectrum as follows:

where $ v_n $ is the flow coefficient and $ \Phi_n $ is the corresponding orientation angle. One can combine $ v_n $ and $ \Phi_n $ together and define a vector notation:

Here $ \{ ...\} $ represents the average over the particle momentum distribution. The dynamical evolution of the bulk matter and the observed anisotropic flows can be successfully described by relativistic hydrodynamics simulations. This implies that the hot and dense QGP matter produced in the heavy-ion collisions at RHIC and the LHC is a strongly-coupled QGP and behaves like a relativistic fluid^[4-8].

Relativistic hydrodynamics is based on the conservation law of energy and momentum, $ \partial_\mu T^{\mu\nu}\!=\!0 $, where $ T^{\mu\nu} $ is the energy-momentum tensor. One can decompose $ T^{\mu\nu} $ as follows:

where $ e $, $ P $, $ u^{\mu} $ are local energy density, pressure and flow velocity of the fluid, respectively, $ \Delta^{\mu\nu}\!=\! g^{\mu\nu} - u^\mu u^\nu $ is the projector transverse to the flow velocity $ u $, $ \pi^{\mu\nu} $ and $ \Pi $ are shear tensor and bulk pressure whose evolution equations can be obtained with the help of the entropy principle or kinetic theory^[26-32]. For Israel-Steward hydrodynamics, we have the following equations for $ \pi^{\mu\nu} $ and $ \Pi $^[26-27]:

where $ \dot{X}\!=\!DX\!=\!(u\cdot \partial) X $ is the comoving derivative on $ X $, $ \theta\!=\!\partial \cdot u $, $ \sigma^{\mu\nu}\!=\!\frac{1}{2}(\nabla^\mu u^\nu + \nabla^\nu u^\mu) - \frac{1}{3} (\nabla \cdot u) \Delta^{\mu\nu} $, $ \nabla^\mu\!=\!\Delta^{\mu\nu}\partial_\nu $, $ \zeta $ and $ \eta $ are the bulk and shear viscosities, $ \tau_{\Pi} $ and $ \tau_{\pi} $ are the corresponding relaxation times. To solve the hydrodynamic equations, the equation of state (EOS) $ P\!=\!P(e) $, which relates the energy density to the pressure, has to be supplemented. The EOS is usually taken from lattice QCD. The production of particles (particlization) from relativistic hydrodynamics is usually performed via Cooper-Frye description^[33]. The produced hadrons continue to interact with each other until the kinetic freeze-out. The hadronic interaction may be simulated via hadron cascade models, such as UrQMD^[34].

There are two important factors in the initial states responsible for the final state anisotropic flows of the particle momentum distribution. The first is the intrinsic collision geometry which is usually an almond shape for a typical non-central collision. Another is the quantum fluctuations of the initial states, which generate various shapes and their combinations for the produced QGP. The initial spatial anisotropies can generate the anisotropies (flows) for the momentum distribution via hydrodynamic evolution. The intrinsic collision geometry is one of the main driving forces to generate elliptic flow $ v_2 $ and other even-order anisotropic flows. The initial state fluctuations not only contributes to even-order harmonic flows, but most importantly can generate triangular flow $ v_3 $ and other odd-order anisotropic flows^[18-22]. The elliptic flow $ v_2 $, triangular flow $ v_3 $ and other higher-order harmonic flows can be well understood by initial state fluctuations and relativistic hydrodynamics, which can not only describe the event-averaged value of $ v_n $, but also the event-by-event distributions of flow $ v_n $ and the correlations between flows of different orders.

Various studies have suggested that the event-by-event fluctuations of flow $ {{V}}_n $ ($ v_n $) can be approximately described by two-dimensional Gaussian distribution (Bessel-Gaussian distribution)^[35-38]. Recently, the deviation of $ v_n $ distributions from the Bessel-Gaussian baseline has been studied using the higher-order moments of $ v_n $ and the comparisons of multi-particle $ v_n $'s^[39-41]. It has been shown by hydrodynamic calculations that the event-averaged $ v_2 $ and $ v_3 $ are linearly correlated with the initial state $ \epsilon_2 $ and $ \epsilon_3 $. However, higher-order $ v_n $'s originate from a combination of initial anisotropy $ \epsilon_n $ and the non-linear response of lower-order flow $ v_m (m<n) $. For example, $ v_4 $ has a non-linear relation to $ v_2 $^[42-43]. Another interesting observation is the negative correlation between $ v_3 $ and $ v_2 $^[42-43], which can be understood from the initial correlation between $ \epsilon_2 $ and $ \epsilon_3 $. Not only the flow magnitudes $ v_n $'s of different orders are correlated^[44], their orientations $ \Phi_n $'s are also correlated with each other^[45-49]. Detailed studies on the event-by-event fluctuations of flow $ v_n $ and the correlations of $ {{V}}_n $ with $ {{V}}_m $ are helpful to achieve more precise determination of the transport properties of QGP.

In heavy-ion collisions, the initial density profiles fluctuate not only in the transverse plane, but also in the longitudinal direction. The longitudinal fluctuations can lead to the fluctuations and decorrelations of flow $ v_n $'s along the pseudorapidity direction^[50-61]. To quantify the effect of longitudinal decorrelations, one may define the following correlation function^{[55, 58]}:

where $ \eta_{\rm{ref}} $ is the reference rapidity. Similarly, one may also define the correlation functions $ r_M(\eta)[n,k] $ and $ r_\Phi(\eta)[n,k] $ using flow magnitude $ v_n $ and orientation $ \hat{\rm{V}}_n $^[59]. Since the above correlation function is, to the leading order, a linear function of $ \eta $, one may simply use its slope to quantify the decorrelation effect. According to the expression $ r[n,k](\eta) \approx 1 - 2 f[n,k] \eta $, one can define the slope parameter as follows:

Similarly, one may also define the slope parameters $ f_M[n,k] $ and $ f_\Phi[n,k] $ for $ r_M(\eta)[n,k] $ and $ r_\Phi(\eta)[n,k] $.

In Fig. 1^[61], we show the effect of longitudinal decorrelation (in terms of the slope parameter $ f $, $ f_M $ and $ f_\Phi $) for Pb+Pb collisions at 5.02 $ A $TeV at the LHC (upper panel) and Au+Au collisions at 200 $ A $GeV at RHIC (middle panel) from a (3+1)-dimensional hydrodynamics calculation using the AMPT^[62] initial conditions. It is interesting to note that within the full AMPT model, significant amount of elliptic flow may come from the escape mechanism^[63]. From Fig. 1, one can see that the longitudinal decorrelation effect is larger at RHIC than at the LHC (see Ref. [64] for STAR result), which may be traced back to the longitudinal structures of the initial states. Since RHIC has lower collision energy and thus has a larger violation of longitudinal boost invariance, the lengths of the initial string structures (in the AMPT model) is smaller than those at the LHC. In the lower panel of Fig. 1, we show our prediction for longitudinal decorrelations in Xe-Xe collisions at 5.44 $ A $TeV at the LHC^[65], which agrees with recent ATLAS data^[66] quite well.

Figure 1.  (color online) The slope parameters $f[n, k]$, $f_M[n, k]$ and $f_\Phi[n, k]$ (with $n\!=\!2,3,4$ and $k\!=\!1$) of the longitudinal decorrelation functions as a function of centrality in Pb+Pb collisions at 5.02 $A{\rm TeV}$ (upper panel)^[61], Au+Au collisions at 200 $A{\rm GeV}$ (middle panel)^[61], and Xe+Xe collisions at 5.44 $A{\rm TeV}$ (lower panel)^[65]. The ATLAS data^{[58, 66]} are compared.

It is worth noting that machine learning and Bayesian analysis have recently been applied to the study of relativistic heavy-ion collisions. Ref. [67] builds a deep convolutional neural network to identify QCD equation of state employed in relativistic hydrodynamics simulations using the simulated final-state particle spectra. Ref. [68] applies Bayesian analysis to a dynamical collision model and a wide variety of experimental data, and simultaneously extracts the estimates of the temperature-dependent specific shear and bulk viscosities. It is found that the specific viscosity of QGP is at least an order of magnitude smaller than that of most common fluids, which implies that QGP behaves like a nearly perfect fluid^[4-8].

Jets that are produced in early hard scatterings have been regarded as valuable probes of QGP in relativistic heavy-ion collisions. During their passage through QGP, jets interact with QGP via elastic and inelastic interactions, and lose some energy during the processes. Such phenomenon is usually called jet quenching^[69-70]. One consequence of jet energy loss is the suppression of large $ p^{}_{\rm{T}} $ hadrons and jets in nucleus-nucleus collisions as compared to the expectation from independent nucleon-nucleon collisions. To quantify such jet quenching effect, one may define the nuclear modification factor $ R_{\rm{AA}} $ as follows:

where $ N_{\rm{coll}} $ is the number of binary nucleon-nucleon collisions. Various experimental results have shown that $ R_{\rm{AA}} $ for high $ p^{}_{\rm{T}} $ hadrons and jets are significantly below unity, which means that jet quenching is mainly due to parton energy loss caused by the final-state interaction between hard jets and hot QGP medium.

To study parton energy loss and jet quenching in relativistic heavy-ion collisions, various theoretical formalisms have been developed for calculating radiative energy loss^[71-76] and collisional energy loss^[77-81] experienced by jet partons inside the dense nuclear matter. Recently, many studies have been performed to improve various approximations used in previous approaches. For example, Ref. [82] extends the Arnold-Moore-Yaffe formalism to include the finite size effect on radiative energy loss. Ref. [83] studies the next-to-eikonal corrections within the path integral formalism. Refs. [84-86] extend the higher-twist formalism by going beyond collinear rescattering expansion and soft gluon emission approximation. Refs. [87-88] reinvestigate the Gyulassy-Levai-Vitev formalism by relaxing soft gluon emission approximation. Refs. [89-90] calculate the correction to multiple soft scattering limit in the Baier-Dokshitzer-Mueller-Peigne-Schiff-Zakharov formalism. In addition, the interference between sequential gluon emissions has also been investigated in Refs. [91-92].

The general framework for the phenomenological studies of jet quenching is based on perturbative QCD factorization theorem, which states that processes involving large momentum transfer can be factorized as the convolution of parton distribution functions (PDFs), hard partonic scattering process, and fragmentation function (FF). For single inclusive high-$ p^{}_{\rm{T}} $ hadron production in proton-proton collisions, the cross section can be expressed as follows^[93-94]:

where $ f(x_a) $, $ f(x_b) $ are PDFs, $ {\rm d}\sigma_{ab \to c} $ is the parton scattering cross section, and $ D_{h/c}(z_c) $ is FF. To study single inclusive high-$ p^{}_{\rm{T}} $ hadron production in heavy-ion collisions, we need to take into account two additional effects. The first is the cold nuclear matter effect due to the modification of PDFs in nuclei. Another is the hot nuclear matter effect caused by the final state jet-medium interaction. Here we provide some details on a linearized Boltzmann transport (LBT) model^[95-101], a quite general framework for simulating heavy and light flavor jet quenching in heavy-ion collisions.

In LBT model, one essentially solves the following Boltzmann equation using Monte-Carlo method:

where $ C_{\rm{el}} $ and $ C_{\rm{inel}} $ represent the collision integrals for elastic and inelastic processes. The elastic scattering between jet partons and medium is simulated using the differential scattering rate calculated from the leading order perturbative QCD matrix elements. The inelastic radiative process for jet partons is simulated using the higher-twist energy loss formalism^{[76, 102-103]}.

Fig. 2 shows the numerical results on nuclear modification factor $ R_{\rm{AA}} $ from Ref. [101], which combines the LBT model with relativistic hydrodynamics to study the energy loss and nuclear modifications of heavy and light flavor jets in the LHC heavy-ion collisions. A next-to-leading order perturbative QCD framework is utilized in order to consistently take into account the contributions from both quarks and gluon to light and heavy flavor hadron productions in high-energy nuclear collisions^[93-94]. The CTEQ5 parameterization^[104] is used for PDFs and the EPS09 parameterization^[105] is used for nuclear shadowing effect. The FFs of charged hadrons, $ D $ meson and $ B $ mesons are taken from Refs. [106-108], respectively. We can see that by including both quark and gluon contributions to light and heavy flavor hadron productions, Ref. [101] obtains a satisfactory description of the $ R_{\rm AA} $ data^[109-112] for light hadrons, $ D $ mesons, $ B $ mesons and $ B $-decayed $ D $ mesons simultaneously over a wide range of $ p^{}_{\rm{T}} $ (8~300 GeV). Below 30~40 GeV, $ B $ meson is less suppressed than charged hadrons and $ D $ mesons due to the large masses of $ b $ quarks. Above 30~40 GeV, Ref. [101] predicts that charged hadrons, $ D $ mesons and $ B $ mesons all have similar $ R_{\rm{AA}} $'s, which may be tested in the future measurements.

Figure 2.  (color online) Nuclear modification factor $R_{\rm{AA}}$ for charged hadrons (a) and $D$ mesons (b) in central 0~10% and 0~80% Pb+Pb collisions at 5.02 $A{\rm TeV}$. Also shown are $R_{\rm{AA}}$ for hadrons initiated by quarks and gluons^[101]. (c) Nuclear modification factor $R_{\rm{AA}}$ for charged hadrons, $D$ mesons, $B$ mesons and $B$-decayed $D$ mesons in central 0~80% (or minimum bias) Pb+Pb collisions at 5.02 $A{\rm TeV}$^[101]. The CMS^[109-112] and ALICE^[113-114] data are compared.

While high-$ p^{}_{\rm{T}} $ heavy hadrons can be utilized to probe the flavor dependence of jet quenching^{[98, 101, 115-121]}, low and intermediate $ p^{}_{\rm{T}} $ heavy hadrons are ideal for studying the hadronization of heavy quarks in heavy-ion environment. Various studies have demonstrated that the coalescence mechanism, in which heavy quarks combine with the thermal partons from the QGP medium to form hadrons, plays important roles for heavy flavor hadron production at low and intermediate $ p^{}_{\rm{T}} $ regime^[122-127]. In particular, coalescence is essential to understand the non-monotonic $ p^{}_{\rm{T}} $ dependence of $ D $ meson $ R_{\rm{AA}} $ from low to intermediate $ p^{}_{\rm{T}} $^[128-129]. Recently, experiments have observed a strong enhancement for $ \varLambda_c $ and $ D_s $ meson yields as compared to that of $ D^0 $ mesons in heavy-ion collisions^[130-132]. Various studies have been performed to explain such enhancement effect^[133-138].

In Ref. [138], a comprehensive coalescence model has been developed including a complete set of both $ s $ and $ p $-wave hadronic states. Such advanced hadronization scheme is then combined with a Langevin-hydrodynamics model^[126-127], which describes elastic and inelastic energy loss of heavy quarks inside the realistic QGP medium. In Figs. 3(a, b), we show the numerical results for $ \varLambda_c/D^0 $ and $ D_s/D^0 $ ratios from Ref. [138]. One can see that the $ p^{}_{\rm{T}} $ differential $ \varLambda_c/D^0 $ and $ D_s/D^0 $ ratios can be described quite well by the coalescence model in Ref. [138]. It is also found that the radial flow of the QGP is essential to understand the chemical composition of the charmed hadrons observed in the experiments. In Fig. 3(c), we show our prediction for the $ B_s/B^+ $ ratio in Pb+Pb collisions at 5.02 $ A $TeV^[139], which agrees with recent CMS data^[140] very well.

Figure 3.  (color online) (a, b) $\varLambda_c/D^0$ and $D_s/D^0$ ratios as a function of $p^{}_{\rm{T}}$ in Au+Au collisions at 200 $A{\rm GeV}$ at RHIC. The results with only fragmentation contribution are also shown^[138]. (c) $B_s/B^+$ ratio as a function of $p^{}_{\rm{T}}$ in Pb+Pb collisions at 2.76 $A{\rm TeV}$ at the LHC^[139]. The STAR data^[131-132] and CMS data^[140] are compared.

In the last decade, fully reconstructed jets have become a hot topic in jet quenching studies in heavy-ion collisions^[141-150]. For full jets, one needs to consider not only the in-medium evolution and energy loss of the leading parton, but also how the shower (radiated) partons interact with the medium and where they go in the final state. For full jet energy loss, one may decompose it into the following components^{[144, 150-152]}:

First, only the out-of-cone radiation is considered as full jet energy loss. Second, some energy inside jet cone may be scattered out of jet cone by elastic interactions. Third, jet partons can deposit energy into the medium or soft shower partons of the jet may thermalize into the medium. Note that some deposited or thermalized energy flows to large angles away from the jets, and some energy is still inside the jet cone. To study the internal structure of full jets, one of the most interesting observables is jet shape function, which describes the energy distribution in the radial $ r $ direction,

where $ R $ is jet radius and $ r \!=\! \sqrt{(\phi \!-\! \Phi_J)^2 \!+\! (\eta - \eta_J)^2} $, with $ \Phi_J $ and $ \eta_J $ being the azimuthal angle and the pseudorapidity of the jet.

In Refs. [151-153], a coupled jet-fluid model has been developed to study the interaction of full jet shower with the dynamically-evolving QGP medium and the nuclear modification of full jet observables in heavy-ion collisions. In this calculation, the evolution of full jet shower is obtained by solving the following Boltzmann equation for the three-dimensional momentum distribution of partons inside the jets $ f(\omega, k_\perp^2, t)\!=\!{\rm{d}}N(\omega, k_\perp^2,t )/{\rm{d}}\omega {\rm{d}}k_\perp^2 $ as follows:

In the above equation, the first and second terms describe the elastic energy loss and transverse momentum broadening effects, respectively, with $ \hat{e} $ and $ \hat{q} $ being the elastic energy loss rate and transverse momentum broadening rate. The third term describes the contribution from medium-induced radiation (splitting). To simulate the medium response to jet-deposited energy and momentum, the hydrodynamic equation with the source term is solved:

where the source term $ J^\nu(x) $ may be constructed from the solution of Boltzmann equation.

Fig. 4 shows the numerical results for full jet shape function and its nuclear modification factor from Refs. [152-153]. One can see that jet-induced medium excitation dominates jet shape function at large $ r $, which means that a lot of jet-deposited energy flows to large angles away from the jet direction. Similar result has also been obtained by other calculations^{[100, 154-156]}. Another interesting result is that for both inclusive jets and $ \gamma $-tagged jets, the nuclear modification of jet shape function for jet $ p^{}_{\rm{T}} > $120 GeV is a non-monotonic function of $ r $, i.e., an enhancement at larger $ r $ and a dip at $ r\sim (0.1\thicksim0.15) $, while the modification factor for $ p^{}_{\rm{T}} > $30 GeV is a monotonic function of $ r $, i.e., an enhancement at larger $ r $ and a dip at $ r \sim 0 $. This implies that the nuclear modification of $ \rho(r) $ has a weak dependence on jet flavors (quark or gluon), but has a strong dependence on jet energies. Such finding may be tested by future measurements with lower energy inclusive jets or high energy $ \gamma $-tagged jets.

Figure 4.  (color online) (a) jet shape function $\rho(r)$ in Pb+Pb collisions at 2.76 $A{\rm TeV}$ at the LHC. The contributions from jet shower and hydrodynamic response are also shown^[152]. (b, c) nuclear modification factor of $\rho(r)$ for single inclusive jets and $\gamma$-tagged jets with two different jet energies in Pb+Pb collisions at 5.02 $A{\rm TeV}$ at the LHC^[153]. The CMS data^[157-159] are compared.

Jet fragmentation function is another interesting observable, describing the distribution of the fractional energy carried by the hadrons inside the jets^[160-161]. Recently, extensive studies have been performed on jet substructure in relativistic heavy-ion collisions. For example, by removing the soft components of the jets using jet grooming algorithms^[162-164], one may reconstruct two hard branches that satisfy certain conditions. One interesting observable is the momentum sharing $ z_{\rm{g}} $ distribution of two subjets, which provides a unique opportunity to study the first hard splitting of the partonic jet in hot and dense QGP and directly probe the medium-induced jet splitting function^[165-175]. In Fig. 5, we show our calculation of the nuclear modification factor for the groomed jet momentum sharing distribution $ p(z_{\rm{g}}) $ in Pb+Pb collisions at 5.02 $ A{\rm TeV} $^[176] for two different jet radius sizes ($ R $=0.2 and $ 0.4 $) and for two different $ z_{\rm{cut}} $ values: $ z_{\rm{cut}} $=0.1 (a) and $ z_{\rm{cut}} $=0.2 (b) compared to recent ALICE data^[177].

Figure 5.  (color online) Nuclear modification factor for groomed jet momentum sharing $z_{\rm{g}}$ distribution in Pb+Pb collisions at 5.02 $A{\rm TeV}$ for two $z_{\rm{cut}}$ values: $z_{\rm{cut}}$=0.1 (a) and $z_{\rm{cut}}$=0.2 (b)^[176]. The ALICE data^[177] are compared.

Small collision systems such as proton-proton and proton-nucleus collisions have been regarded as the references for studying various signatures and properties of QGP in nucleus-nucleus collisions. Intriguingly, strong collectivity has been observed in proton-lead^[178-180], deuteron-gold^[181], and even high multiplicity proton-proton collisions^[182]. The origin of the collectivity in small collision systems has now become one of the most interesting and important topics in heavy-ion physics. The final answer to it may significantly affect our understanding of the collective behaviors and transport properties of the dense QCD matter produced in nucleus-nucleus collision systems.

One of the most popular explanations of the collectivity in small collision systems is the hydrodynamics approach, in which the anisotropic flows are viewed as a final state effect generated from relativistic hydrodynamic evolution of the fireball with initial state fluctuations and anisotropies^[183-189]. It is usually assumed by hydrodynamics approach that the QGP fluid is formed. However, no jet quenching has been observed in small collision systems so far, which may be understood since the size of the QGP fireball is small. Another alternative explanation is based on the color glass condensate (CGC) framework or the saturation dynamics which shows that the correlations between the partons originated from the projectile proton and the dense gluons inside the target nucleus can also provide significant collectivity for light hadrons^[190-193].

Recently, experiments at the LHC have also observed sizable collectivity for charm quarks in small collision systems: large elliptic flow $ v_2 $ for both $ D^0 $ and $ J/\Psi $ mesons in proton-lead collisions^[194-196], though the $ v_2 $ values are smaller than $ v_2 $ of light hadrons. While the final state effect from hydrodynamics evolution can provide successful description for light hadron $ v_2 $ in small collision systems, it is usually difficult for hydrodynamics to generate large collectivity for heavy mesons. This is simply because: heavy quarks with large masses in general do not flow as much as light quarks and gluons. As shown by Ref. [197], the final state effects caused by the formation of an expanding QGP medium cannot account for the large $ v_2 $ values for $ J/\Psi $ observed in the LHC proton-lead collisions.

In Ref. [198], the elliptic flow $ v_2 $ of $ J/\Psi $ is studied within a simplified model based on the color evaporation model and the dilute-dense factorization in the CGC framework. It considers a gluon and a quark originated from the projectile proton, with the quark serving as the reference whereas the gluon splitting into a pair of heavy quarks. During their passage through the target nucleus, the incoming partons interact with the dense gluons in the nucleus, which generates correlations between the heavy quark pair and the accompanying reference quark. As a result, a significant amount of elliptic flow $ v_2 $ for $ J/\Psi $ can be generated due to the correlations between $ J/\Psi $ and light quarks (hadrons), as shown in Fig. 6. This calculation also predicts that heavier $ \Upsilon $ has similar $ v_2 $ values as $ J/\Psi $, which may be tested in the future measurements.

Figure 6.  (color online) The integrated $v_2$ of $J/\Psi$ and $\Upsilon$ compared to the reference light quark $v_2$ as a function of the saturation scale $Q_{\rm s}^2$ (a) and the transverse momentum dependence of $v_2$ of $J/\Psi$ and $\Upsilon$ (b) in p+Pb collisions at $\sqrt{s^{}_{\rm{NN}}}$ = 8.16 TeV at the LHC^[198]. The CMS data^[199] are compared.

One of the main purposes of relativistic heavy-ion collisions is to create the strong-interaction matter at various temperatures and densities, and to study the novel properties of quark-gluon plasma. In the last decades, tremendous progresses have been made in the field of high-energy nuclear collisions. In this paper, we have reviewed some of the recent developments in the soft and hard probes of QGP related to the heavy-ion experiments at RHIC and the LHC. Our main focus was given to anisotropic collective flows and jet quenching, which have been regarded as two most important signatures of QGP in relativistic heavy-ion collisions.

For collective behaviors of the QCD matter, we have discussed initial-state fluctuations, hydrodynamical evolution and final-state flows, correlations and longitudinal decorrelations. One important direction is to utilize the advanced statistical tools to analyze various bulk observables for various collision systems at various energies, which can help us to obtain more precise determination of the transport properties, such as the viscosities, of the QGP. Regarding hard probes of the QGP, we have discussed the energy and flavor dependences of jet quenching, heavy quark hadronization, full jet structure/substructure and medium response. With more realistic simulation of jet-medium interaction and more detailed analysis of jet-related observables, one should be able to probe the properties and internal structure of the QCD matter at various energy and length scales. For collectivity in small collision systems, we have discussed the observed collective flows of light and heavy flavor hadrons in proton-nucleus collisions. One future task is to disentangle the contribution from initial state correlations and that from the final state interactions, which is also essential for understanding the origins of anisotropic flows in large nucleus-nucleus collision systems. The ultimate goal is to map out the novel properties and phase structures of the QCD matter at various conditions and achieve a comprehensive understanding on the non-equilibrium dynamics of high-energy nuclear collisions.