**1. Introduction**

The study of particle production in relativistic heavy-ion collisions is important for understanding the dynamics of the evolution of hot and dense systems formed in these processes. The investigation of such complicated phenomena is much more effective when the analysis of the experimental data is accompanied by simulations within a realistic model of the collision. Theoretical analysis often allows finding the most adequate interpretation of the measured data, clarifying the role of specific factors in the formation of observed results, and understanding the interplay (often quite complicated) between these factors.

In the recent papers [1–3], the integrated hydrokinetic model (iHKM) [4] was applied to the description of particle production in Au+Au collisions at the top RHIC energy √*sNN* = 200 GeV and Pb+Pb collisions at the LHC energies √*sNN* = 2.76 TeV and √*sNN* = 5.02 TeV. The analysis of the particle momentum spectra, particle yields, and particle number ratios shows in particular that a successful description of the data can be reached in the model even using different equations of state for quark-gluon matter at the hydrodynamics stage of expansion, if the initial energy-density profile is correspondingly rescaled. The particle number ratios, calculated with inelastic processes switched off at the final "afterburner" stage of the collision, differ from the same ratios calculated in "full mode" (i.e., when the inelastic processes are taken into account) and, unlike the latter, depend on the utilized equation of state. This fact shows the importance of the afterburner particle scatterings for the formation of the observed ratio values.

In the article [5], a conclusion about grea<sup>t</sup> influence of intensive particle interactions with the hadronic medium at the late stage of the collision on the number of detected *K*∗(892) resonances (which were reconstructed in the experiment through the products of their decay *K*∗(892) → *Kπ*) was made based on the results of iHKM simulations. The numerous hadronic interactions, leading to multiple rescatterings and recombinations of resonance decay products, also resulted in larger estimated times of maximal emission (Since in iHKM, we do not have a sudden kinetic freeze-out, but particles emit continuously from the system, we utilize the time of maximal emission *τm*.*e*. concept in our studies. For particles of a certain species and from a certain *pT* bin, *τm*.*e*. defines a spacelike hypersurface *τ* = *τm*.*e*. with a thin 4D layer, adjacent to it, from which most of considered particles are emitted.) for kaons as compared to pions, as showed the analysis of the particle emission pictures, obtained from iHKM in the considered cases [1,6].

The mentioned results sugges<sup>t</sup> that chemical and kinetic "freeze-out" in relativistic A+A collisions is rather continuous, not sudden. Such a conclusion looks reasonable, since sudden chemical freeze-out would mean an abrupt switching from the expansion of chemically equilibrated matter with very intensive inelastic reactions to the evolution of the system, where inelastic reactions are absent. As for a sudden kinetic freeze-out, it would suppose a sudden switching from very large hadron-hadron interaction cross-sections (in a nearly perfect hydrodynamics regime) to zero cross-sections (specific for the free streaming particles case). Both sharp transitions can hardly be expected from the theoretical point of view.

In the present work, we aim to apply iHKM to the description of particle production in Xe+Xe collisions at the LHC energy √*sNN* = 5.44 TeV and find out if the results will be similar to those obtained in our previous studies.

#### **2. Results and Discussion**

In order to simulate the process of matter evolution during a high-energy heavy-ion collision, we utilized the integrated hydrokinetic model [2,4]. It consists of several blocks, each modeling a certain stage of collision. The first block deals with the early prethermal system's dynamics, in the course of which it gradually transforms from the initial, non-thermalized state to the state of local equilibrium (chemical and thermal). To describe this process in iHKM, we simulated the evolution of the energy-momentum tensor of the system within an energy-momentum transport approach in a relaxation time approximation.

To obtain the initial conditions (i.e., the initial non-thermal energy-momentum tensor) for the pre-equilibrium stage, we used a combined approach: our tensor *<sup>T</sup>μν*(*x*) is defined by the initial distribution function *f*(*<sup>x</sup>*, *p*), so that *<sup>T</sup>μν*(*x*) = *d*3 *p f*(*<sup>x</sup>*, *p*)*p<sup>μ</sup> p<sup>ν</sup>*/*p*0. The *f*(*<sup>x</sup>*, *p*) can be factorized into space-time and momentum parts. The form of the momentum part of the function is based on the color glass condensate approach:

$$f\_0(p) = g \exp\left(-\sqrt{\frac{(p \cdot \mathcal{U})^2 - (p \cdot \mathcal{V})^2}{\lambda\_\perp^2} + \frac{(p \cdot \mathcal{V})^2}{\lambda\_\parallel^2}}\right),\tag{1}$$

where *Uμ* = (cosh *η*, 0, 0, sinh *η*), *Vμ* = (sinh *η*, 0, 0, cosh *η*), *η* is space-time rapidity, and initial momentum anisotropy Λ = *<sup>λ</sup>*⊥/*λ* = 100 [4]. The spatial part of the initial distribution function is generated by the GLISSANDOcode [7], which implements the Monte Carlo Glauber model. It gives us the initial transverse energy-density profiles (*b*,**r***T*) (here, **r***T* is the transverse plane radius vector and *b* is the impact parameter, denoting the collision centrality class):

$$\epsilon(b, \mathbf{r}\_T) = \epsilon\_0(\tau\_0) \frac{(1 - a)N\_{\text{w}}(b, \mathbf{r}\_T)/2 + aN\_{\text{bin}}(b, \mathbf{r}\_T)}{(1 - a)N\_{\text{w}}(b = 0, \mathbf{r}\_T = 0)/2 + aN\_{\text{bin}}(b = 0, \mathbf{r}\_T = 0)}. \tag{2}$$

The coefficient 0(*<sup>τ</sup>*0) here is the main iHKM parameter, which allows adjusting the model to the description of different types of collisions, the initial energy density in the center of the system at the initial proper time being *τ*0. The other parameter *α* defines the proportion between the contributions *Nbin*(*b*,**r***T*) and *Nw*(*b*,**r***T*) from the "binary collisions" and the "wounded nucleons" models to the GLISSANDO energy-density profile. Note, that the quadrupole deformation of somewhat prolate xenon nuclei with deformation parameter *β*2 = (0.18 ± 0.02) [8] also can be taken into account within the GLISSANDO code (for this purpose, the spherical nuclear radius *R*0 is substituted in simulations

by the polar angle *θ*-dependent value *R*(*θ*) = *<sup>R</sup>*0[<sup>1</sup> + *β*2*Y*<sup>0</sup> 2 (*θ*) + *β*4*Y*<sup>0</sup> 4 (*θ*)]), where *Y*<sup>0</sup> 2 , *Y*<sup>0</sup> 4 are spherical harmonics. For the current study, we used the GLISSANDO profiles, generated with *β*2 = 0.18.

After the prethermal stage follows the second iHKM block, which describes the hydrodynamical expansion of locally thermalized matter, using the relativistic viscous hydrodynamics approach based on the Israel–Stewart formalism. The energy-momentum tensor *<sup>T</sup>μν*(*x*) of the system now acquires the hydrodynamical form and gives us the hydrodynamical values of local energy density (*x*), pressure *p*(*x*), and collective velocity *<sup>u</sup><sup>μ</sup>*(*x*). The shear viscosity to entropy ratio *η*/*S* is assumed to equal its minimal theoretical value 1/(<sup>4</sup>*π*) ≈ 0.08, which was calibrated in [2] as giving the best description of experimental data for Pb+Pb LHC collisions at 2.76 *A* TeV.

The third stage in iHKM describes the particlization of the system, when we switch from describing it in terms of a continuous medium to the description in terms of particles. We assumed that the particlization took place at the isotherm hypersurface with temperature *Tp*, which in the general case depends on the used equation of state (EoS) for quark-gluon matter. The latter is matched at *T* = *Tp* with the EoS of an ideal chemically equilibrated hadron-resonance gas, consisting of *N* = 329 well-established hadron states.

The switching hypersurface *σsw* is constructed in the course of the hydrodynamic evolution of the system using the Cornelius code [9–11]. The Cooper–Frye prescription is applied to convert the fluid to particles:

$$\left.p^{0}\frac{d^{3}N\_{i}(\mathbf{x})}{d^{3}p}\right|\_{d\sigma(\mathbf{x})} = d\sigma\_{\mu}(\mathbf{x})p^{\mu}f\_{i}(p\cdot\mathbf{u}(\mathbf{x}),T(\mathbf{x}),\mu\_{i}(\mathbf{x})),\tag{3}$$

where *dσ*(*x*) is the particlization hypersurface element near point *x* and *i* enumerates the particle species. The Grad ansatz is utilized to account for the viscous corrections to the equilibrium distribution function *f eq i* (*<sup>x</sup>*, *p*) [12]. The corrections are assumed to be the same for all hadron species. Thus, the Equation (3) can be rewritten (in the fluid local rest frame) as:

$$\frac{d^3 \Delta N\_i}{dp^\* d(\cos \theta) d\phi} = \frac{\Delta \sigma\_{\mu}^\* p^{\*\mu}}{p^{\*0}} p^{\*2} f\_i^{\neq 0} \left( p^{\*0}; T, \mu\_i \right) \left[ 1 + (1 \mp f\_i^{\neq 0}) \frac{p\_\mu^\* p\_\nu^\* \pi^{\*\mu \nu}}{2T^2(\epsilon + p)} \right],\tag{4}$$

where *π*<sup>∗</sup>*μν* is the shear stress tensor.

The particle coordinates and momenta are generated according to distributions (4) using the Monte Carlo procedure, so that the total hadron multiplicity in each event is randomly sampled in accordance with the Poisson distribution with a mean value equal to the mean total number of particles *Ntot* at the switching hypersurface *σsw*, which is calculated based on (4). The sort for each generated particle is chosen randomly with probability *Ni*/*Ntot*, where *Ni* is the mean number of hadrons of sort *i* at *σsw* according to (4).

At the next, final stage of the evolution, the generated particles, which now form our system, undergo multiple elastic and inelastic scatterings, and all possible resonance decays take place. This final block in iHKM was realized with the help of the UrQMDmodel [13,14].

As compared to other collision models utilizing the hybrid (hydro+cascade) approach (see, e.g., [15–19]), which succeed in reproducing the experimental data on soft observables in high-energy heavy-ion collisions, iHKM has an important advantage: the first prethermal stage of the matter evolution in iHKM allows one to describe the development of collective flow and other effects adequately, which starts due to the finiteness and azimuthal asymmetry of the system already at the early initial time *τ*0 ≈ 0.1 fm/*<sup>c</sup>*, when the matter is not ye<sup>t</sup> thermalized. This allows iHKM to describe all main bulk observables together with the femtoscopy radii for different collision types simultaneously.

The model parameter values that allow describing Xe+Xe collisions at 5.44 *A* TeV were fixed as those giving the best agreemen<sup>t</sup> with the experimental mean charged particle multiplicity dependency on centrality (see Figure 1) and the pion spectrum for the most central events. The 0 = 445 GeV/fm<sup>3</sup> at *τ*0 = 0.1 fm/*c* and *α* = 0.44 values were obtained from these fits. The chosen maximal energy density value 0 was between that for the top RHIC energy Au+Au collisions (235 GeV/fm3) and the one for 2.76*A* TeV Pb+Pb LHC collisions (679 GeV/fm3). The Xe+Xe value of *α* was the maximal among those used in our studies so far (for top RHIC energy case *α* = 0.18 and for both Pb+Pb LHC cases *α* = 0.24). A relatively large *α* value for Xe+Xe collisions as compared to the Pb+Pb collisions with close energy, which corresponds to a larger binary collision contribution to the initial energy-density distribution, can be connected with a smaller size of the Xe nucleus in comparison with the Pb nucleus. For quark-gluon phase at the hydrodynamics stage, the HotQCD Collaboration's equation of state was applied in current analysis [20] with particlization temperature *Tp* = 156 MeV.

**Figure 1.** The mean charged particle multiplicity *dN*ch/*dη* for |*η*| < 0.5 in Xe+Xe collisions at the LHC energy √*sNN* = 5.44 TeV for different centrality classes calculated in the integrated hydrokinetic model (iHKM) and measured by the ALICE Collaboration [8].

In Figures 2–5, the transverse momentum spectra for all charged particles, pions, kaons, and protons calculated in the integrated hydrokinetic model for eight centrality classes are demonstrated together with the corresponding experimental data from the ALICE Collaboration [8,21]. From the plots, one can see that the model describes the experimental *pT* spectra for all the particle species and all the centrality classes quite well. Note in particular that iHKM reproduces the measured pion spectra in a wide *pT* interval, including the region of soft momenta. A similar description quality was reached earlier in [3,4] for 2.76*A* TeV and 5.02*A* TeV Pb+Pb collisions at the LHC.

**Figure 2.** The iHKM results on transverse momentum spectra of all charged particles compared to the experimental data from the ALICE Collaboration [8] for Xe+Xe collisions at the LHC energy √*sNN* = 5.44 TeV. The datasets for different centralities are scaled for better visibility, |*η*| < 0.8.

 **Figure 3.** The same as in Figure 2 for pions, |*y*| < 0.5 (ALICE Collaboration points are taken from [21]).

 **Figure 4.** The same as in Figure 2 for kaons, |*y*| < 0.5 (ALICE Collaboration points are taken from [21]).

**Figure 5.** The same as in Figure 2 for protons, |*y*| < 0.5 (ALICE Collaboration points are taken from [21]).

In Figure 6, one can see the results of the calculation in iHKM for different particle number ratios in comparison with the experimental data from the ALICE Collaboration [22–24]. The points correspond to central collision events of Xe+Xe collisions at the LHC energy 5.44*A* TeV and previously analyzed within iHKM Pb+Pb collisions at the LHC energies 2.76*A* TeV and 5.02*A* TeV. The shown model points, corresponding to full iHKM calculation, which included inelastic processes at the afterburner stage, were in good agreemen<sup>t</sup> with the data. The points for reduced regime without inelastic processes in most cases gave a worse or totally inadequate description of the experiment and thus are not shown. Most ratios of given particle yields for different collision types were close to each other, except for the *φ*/*K* ratio for the LHC energy 2.76 *A* TeV. A possible reason for a lower value of the *φ*/*K* ratio in 2.76 *A* TeV collisions could be the shorter and less intensive "afterburner" stage of matter evolution in this case, at which multiple particle scatterings took place, leading in particular to recombination of *φ* resonances from their decay products (kaons). This effect led to an increase of the observed *φ* yields, and so, if in the 5.02 *A* TeV and 5.44 *A* TeV cases, it is more pronounced, the corresponding *φ*/*K* ratios will be larger.

In general, summarizing the presented results, we concluded that iHKM successfully described the particle production in Xe+Xe collisions at the LHC energy √*sNN* = 5.44 TeV, similarly to Pb+Pb collisions at the LHC and Au+Au collisions for RHIC cases, described in previous works. The final stage of the collision played an important role in the formation of bulk observables and supported the suggestion about the continuous character of the chemical and kinetic freeze-out.

**Figure 6.** The iHKM results on various particle number ratios for central Pb+Pb collisions at the LHC energies √*sNN* = 2.76 TeV and √*sNN* = 5.02 TeV and Xe+Xe collisions at the LHC energy √*sNN* = 5.44 TeV. The model points are compared to the experimental ones, where measured data are available [22–24].

**Author Contributions:** Conceptualization, Y.S.; data curation, M.A. and V.S.; funding acquisition, Y.S.; investigation, M.A. and V.S.; methodology, Y.S.; project administration, Y.S.; supervision, Y.S.; visualization, M.A. and V.S.; writing, original draft, V.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was carried out within the scope of the International Research Network "EUREA: European Ultra Relativistic Energies Agreement," and the corresponding Agreement F20-2020 with the National Academy of Sciences (NAS) of Ukraine. The work was partially supported by the Tomsk State University Competitiveness Improvement Program.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.
