Centrality Selection Effect on Elliptic Flow Measurements in Relativistic Heavy-Ion Collisions at NICA Energies

: The elliptic ﬂow ( v 2 ) of produced particles is one of the important observables sensitive to the transport properties of the strongly interacting matter created in relativistic heavy-ion collisions. Detailed differential measurements of v 2 are also foreseen in the future Multi-Purpose Detector (MPD) experiment at the Nuclotron based Ion Collider fAcility (NICA) at collision energies √ s NN = 4–11 GeV. Elliptic ﬂow strongly depends on collision geometry, deﬁned by the impact parameter b . Usually b is an input to theoretical calculations and can be deduced from experimental observables in the ﬁnal state using the centrality procedure. In this work, we investigate the inﬂuence of the choice of centrality procedure on the elliptic ﬂow measurements at NICA energies.


Introduction
The main purpose of the MPD detector at NICA is to explore the phase diagram of the strongly interacting matter in the region of high net-baryon densities [1]. Anisotropic collective flow is one of the most important observables to probe the equation of state and transport properties of matter created in relativistic heavy-ion collisions [2]. It can be quantified according to the Fourier coefficients v n in the expansion of the particle azimuthal distribution relative to the collision symmetry plane given by the angle Ψ n : dN/dφ ∝ 1 + ∑ n=1 2v n cos(n(φ − Ψ n )), (1) where n is the order of the harmonic, and φ is the azimuthal angle of a particle. The coefficients v n can be calculated as v n = cos[n(ϕ − Ψ n )] , where the brackets denote the average over the particles and events. The elliptic flow (v 2 ) is one of the most extensively studied observables in relativistic nucleus-nucleus collisions and it has been measured in different experiments in the last three decades [3]. The detailed measurements of v 2 of various hadrons produced in relativistic heavy-ion collisions at the Relativistic Heavy-Ion Collider (RHIC) and at the Large Hadron Collider (LHC) played an essential role in the discovery of the strongly coupled Quark-Gluon Matter (QGM) [4,5]. The elliptic flow signal strongly depends on collision geometry and viscous hydrodynamic studies suggest that v 2 stems from the evolution of the medium in the presence of initial-state anisotropies, determined by the eccentricity ε 2 of the overlap zone in nucleus-nucleus collisions. The v 2 signal is recognized to be almost linearly correlated to ε 2 : v 2 = k 2 × ε 2 [5][6][7]. The coefficient k 2 encodes the medium response, which is sensitive to the temperature (T) dependent specific shear viscosity (i.e., the ratio of shear viscosity to entropy density η/s(T)) of strongly interacting system produced in the collisions [6][7][8]. Thus, the comparison of viscous hydrodynamical model calculations to the v 2 /ε 2 for various collision centralities is commonly employed to estimate the average η/s(T) over the system's evolution [5,8]. The impact parameter b, defined as the transverse distance between the centers of the two colliding nuclei, is the best variable that can be used to define the collision centrality and a key input to most theoretical calculations of heavy-ion collisions. However, it can not be directly measured in experiments [9]. There, the heavy-ion collisions can be characterized by the measured multiplicities N ch of the produced particles or by the forward (backward) rapidity energy E sp , which is an approximation for the number of projectile (target) spectators. The centrality procedure is based on the correlation between measured N ch (E sp ) and b, which can be deduced by fitting a specific model of particle production to the experimental data [10]. Usually, the correlation between the impact parameter b and the multiplicity N ch is determined using the Monte-Carlo Glauber (MC Glauber) method combined with a simple particle production model [11]. The modeled multiplicity is assumed to be a function of the number of participating nucleons (N part ) and the number of binary interactions between nucleons (N coll ), which one obtains from the output of the MC Glauber model. The particle multiplicity distribution N f it ch can then be fitted to the experimentally measured one [12,13]. Centrality classes are defined by sharp cuts on N ch and the corresponding mean values of b for each class determined from MC Glauber events. Although this approach seems to be well established, it may suffer from large systematic uncertainties at low multiplicities and assumptions about the particle production mechanism [10]. Recently, a new modelindependent (Γ-fit) method for reconstructing the impact parameter distributions was proposed [14]. The main assumption is that the fluctuations of the N ch used to determine the centrality at a fixed impact parameter b follow a gamma distribution [15].
In this work, we investigate the influence of the choice of centrality procedure on the elliptic flow measurements at NICA energies: √ s NN = 4-11 GeV [16]. To address this issue, Au+Au collisions were simulated using Ultra-relativistic Quantum Molecular Dynamics (UrQMD) [17] and A Multi-Phase Transport (AMPT) [18] models. The elliptic flow results from events binning with the Γ-fit, and MC Glauber methods are compared to those with the true impact parameter from the models. The paper is organized as follows. Section 2 briefly describes the used transport models (UrQMD and AMPT), while Section 3 discusses the MC Glauber and Γ-fit methods for centrality determination. Section 4 discusses the procedures used to determine the elliptic flow coefficient v 2 . The effects of centrality filters on v 2 results are presented and discussed in Section 5. Finally, a summary is given in Section 6.

A Brief Description of the Models
In this work, the relativistic Au+Au collisions at √ s NN = 5, 7.7 and 11.5 GeV are simulated by the string melting (SM) version of the AMPT [18] [20][21][22]. Figure 1 shows the comparison of the p T -differential v 2 of charged hadrons between different transport models and STAR published data [20] (blue solid circles) in the 20-30% central Au+Au collisions at √ s NN = 7.7 GeV. The event plane and centrality in the analysis of the model events were determined in the same way as in the real STAR data analysis. An elliptic flow analysis was also performed in the same way using the η-sub event plane method. See [16,20] for further details. The left panel of Figure 1 represents the results for hybrid models with QGM formation: 3D viscous hydro + hadronic cascade vHLLE+UrQMD model [23,24] (open circles) and string melting version of AMPT SM [18] (open boxes). Both models provide a relatively good description of the published v 2 (p T ) results of charged hadrons.
However, the following models: UrQMD [17], DCM-QGSM-SMM [25], JAM [26], and SMASH [27], which only take the hadronic interactions into account (no QGM formation), significantly underestimate the published v 2 values. See the right panel of Figure 1 for more At an energy √ s NN = 4.5 GeV, the hadronic cascade mode of the JAM and UrQMD models can provide a relatively good description of the published STAR v 2 (p T ) results for charged pions and protons. See [19] for the details. This may indicate that at energies √ s NN ≤ 4-6 GeV, the hadron gas phase dominates and models based on only hadron transport can be used to describe the v 2 measurements. Based on these results, we used AMPT SM and the cascade mode of UrQMD to cover the full energy range √ s NN = 4-11 GeV. The AMPT model [18] is a hybrid model, with the fluctuating initial conditions based on the Heavy-Ion Jet Interaction Generator (HIJING) two-component model. In the string melting version of the AMPT SM model, which is used in the present study, hadrons produced from excited strings in the HIJING model are converted to their valence quarks and antiquarks, and the space-time evolution of QGM evolution is then modeled using Zhang's parton cascade (ZPC) model. At hadronization, quarks and antiquarks in the AMPT model are converted to hadrons via a quark coalescence model; following this, the hadronic interactions are modelled by A Relativistic Transport (ART). We have generated events with the AMPT-SM (version v2.26t7) model with the partonic cross section: σ p = 1.5 mb, which qualitatively describes the STAR data of v 2 (p T ) of charged hadrons at √ s NN = 7.7 and 11.5 GeV, Figure 2. The Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model [17] is a microscopic transport approach based on the binary elastic and inelastic scattering of hadrons, resonance excitations and decays as well as string dynamics and strangeness exchange reactions. We used version 3.4 of the UrQMD with the default set of parameters in the cascade mode.
For each model, a Monte Carlo event sample of 50 millions minimum bias Au + Au collisions has been generated for collision energies √ s NN = 5, 7.7, and 11.5 GeV.

Multiplicity-Based Centrality Determination
In this section, we briefly discuss the procedures to determine the centrality of collisions with the Multi-Purpose Detector (MPD) at NICA. The MPD was designed as a 4π spectrometer for detecting charged hadrons, electrons and photons in heavy-ion collisions at high luminosity. In the first stage of operation in 2024, the MPD will consist of the Time Projection Chamber (TPC), the Time-Of-Flight (TOF) detector, the electromagnetic calorimeter (ECal), and the forward hadron calorimeter (FHCal). See the left panel of Figure 3 for further information. TPC will provide 3D tracking of charged particles, as well as measuring the specific ionization energy loss dE/dx to identify the particles with |η| < 1.2. More details about the detector subsystems of MPD and their performance can be found in [1]. The uncorrected multiplicity of charged particles N ch in the TPC (MPD) is used to determine and define the centrality classes. The pseudorapidity region |η| < 0.5 of TPC (MPD) used for centrality selection is marked by the green vertical band. See the right part of Figure 3 for further details. The |η| range for the multiplicity of charged particles is similar to that used by the STAR experiment for centrality definition during the first Beam Energy Scan program at RHIC: √ s NN = 7.7-62.4 GeV [20,22]. The definition of the centrality classes are based on the application of MC Glauber [12,13,28] and Γ-fit [14,15,28] methods.
The present analysis is based on the implementation of the Monte Carlo version of the Glauber model, as described in refs. [11,12]. An input of the MC Glauber model is the nucleon density ρ(r) inside the nucleus. It is usually parametrized by Fermi distribution: where R is the radius of the nucleus (R = R Au = 6.55 ± 0.05 fm for 197 Au nucleus), the constant ρ 0 corresponds to the density at the center of the nucleus. The skin thickness of the nucleus a defines how abruptly the density falls at the edge of the nucleus (a = 0.544 ± 0.010 fm). The nucleus-nucleus collision is treated as a sequence of independent binary nucleon-nucleon collisions, where the nucleons travel on straightline trajectories and the inelastic nucleon-nucleon cross section σ inel NN is assumed to depend only on the collision energy. Two nucleons from different nuclei are assumed to collide if the relative transverse distance d between centers is less than the distance corresponding to the inelastic nucleon-nucleon cross section: d < σ inel NN /π. For selected energies, the values of σ inel NN are set to 29.4, 29.7, and 31.3 mb for √ s NN = 5, 7.7 and 11.5 GeV, correspondingly [29].
The output of the MC Glauber model includes the geometrical properties of the simulated collisions: the impact parameter b, number of binary nucleon-nucleon collisions (N coll ), and the number of participating nucleons (N part ). The procedure for centrality determination includes fitting experimentally measured particle multiplicity N ch with an MC Glauber [12,13,28] : where P µ,k is the negative binomial distribution (NBD) with mean µ and width k. N a ( f ) is a number of ancestors (number of independent sources), f characterizes the fraction of hard processes, N part and N coll are the number of participants and the number of binary collisions from MC Glauber model output. The optimal set of parameters f , µ and k can be found from the minimization procedure applied to find the minimal value of the χ 2 , which is defined as follows: where F i f it and F i data are values of the fit function and fitted histogram at a given bin i, ∆F i f it and ∆F i data are corresponding uncertainties, and n low and n high are the lowest and highest fitting ranges, correspondingly.
For the fit procedure, a 20 × 10 6 MC Glauber events were generated for each energy point. A grid of k and f parameters was formed with corresponding χ 2 values for each (k, f ) combination: k ∈ [1, 50] with a step of 1 and f ∈ [0, 1] with a step of 0.01.
As an example, Figure 4 shows the charged particle multiplicity N ch distribution (open squares) for Au+Au collisions at √ s NN = 7.7 GeV from the UrQMD (left) and AMPT SM (right) models with overlay fitted distribution using the MC Glauber approach (blue solid triangles). The multiplicities N ch below 15 for UrQMD and 10 for AMPT-SM model events were excluded from the fit. This defines the so-called "anchor point" below which the centrality determination is not reliable. The ratio (N f it ch /N ch ) of the fit to the data shows the quality of the procedure, see the bottom part of Figure 4. After finding the optimal set of the fit parameters, one can easily estimate the total cross-section and all events can be divided into groups with a given range of total cross-section (0-10%, 10-20% etc). See the black solid vertical lines in Figure 4 for more information. High multiplicity events have a low average b (central collisions) and low multiplicity events have a large average b (peripheral collisions). For each centrality class, the mean value of the impact parameter b and its corresponding standard deviation was found using simulated information from MC Glauber model events. Figure 5 shows the centrality dependence of b for UrQMD (left) and AMPT SM (right) model events denoted by open symbols. The b from the MC Glauber approach (closed symbols) are presented for comparison. The average impact parameter and the width of its distribution estimated with the MC Glauber approach are consistent with the values used in the models by 3-4%, see the bottom plots in Figure 5.  In contrast to the MC Glauber method, the Γ-fit method does not require any modeling of the collision dynamics and can be used over a broad range of collision energies: from √ s NN = 5.44 TeV [14,15] to the bombarding energy of 25 AMeV [30]. The method shows that the problem of reconstructing impact parameter b from the measured multiplicity N ch is a typical inverse problem, which can be solved using a deconvolution method. The main element is the fluctuation kernel which is used to model multiplicity fluctuations P(N ch |b) at a fixed impact parameter b. The fluctuations of the multiplicity can be described by the gamma distribution [14,15,28]: where Γ(k) is the gamma function and two parameters k(b) and θ(b) corresponding to the mean, N ch , and to the variance, σ N ch : N ch = kθ, σ N ch = √ kθ. Similar to the multiplicity N ch , which is always positive, the gamma distribution is only defined for N ch ≥ 0. It can be considered as a continuous version of the negative binomial distribution (NBD), which has long been used to fit multiplicity distributions in heavy-ion collisions [14,15,30]. The normalized measured multiplicity distribution, P(N ch ), can be obtained by summing the contributions to multiplicity at all impact parameters: where P(b) is the probability distribution of the impact parameter, and c b denotes the centrality: depends on the probability P inel (b) of an inelastic collision occurring at given b, and σ inel is the inelastic nucleus-nucleus cross section. P inel (b) 1 and c b πb 2 /σ inel , except for peripheral collisions. For the variable k, one can use the following parameterization: we fit P(N ch ) to the experimental distribution of N ch using Equations (5) and (6) [14,15,28]. The fit for the reconstruction of the probability of N ch at fixed c b : P(N ch |c b ).
The fitting procedure has been tested for the same charged particle multiplicity N ch distribution from the UrQMD (left) and AMPT SM (right) models, see Figure 6. The result of the Γ-fit is shown as red solid circles. The value of σ inel = 677 fm 2 was used for the Au+Au collision system. In the fit, we exclude values of N ch < 15 below the "anchor point" where a fraction of events are missed. We normalized P(N ch ) in such a way that the fraction of events above the anchor point matches the measured value of N ch . The bottom plots in Figure 6 show the ratio of the resulting fit functions to the charged particle multiplicity distribution. The ratio plots show that the Γ-fit method can reproduce the charged particle multiplicity distribution with good accuracy. Once the probability of N ch at a fixed c b is reconstructed, the probability distribution of b, at a fixed N ch , can be extracted using Bayes' theorem: P(b|N ch ) = P(N ch |b)P(b)/P(N ch ), where P(N ch |b) = P(N ch |c b ) and c b πb 2 /σ inel [14,15,28]. Extending this reconstruction to a finite centrality bin, corresponding to an interval N low ch < N ch < N high ch is straightforward upon integration over N ch :    Results for the Γ-fit approach tend to be in better agreement with model data. However, it should be noted that this approach requires the total integral of the multiplicity distribution to be evaluated separately. Thus, the Γ-fit method is more sensitive to any bias, such as trigger inefficiencies, that could distort the estimation of the total integral of the multiplicity distribution.

Methods for Elliptic Flow Measurements
A significant part of the published elliptic flow measurements by the STAR experiment at RHIC, including the Au+Au beam energy scan programs, have been performed using the traditional events plane and Q-cumulant methods [20][21][22]. In both cases, the particles detected in the TPC of the STAR experiment (|η| < 1.0) have been used. In this section, we briefly discuss how the event plane and Q-cumulant methods can be used for the measurements of elliptic flow of the produced particles with an MPD detector system at NICA [31]. Similar to the STAR experiment, the particles detected in the TPC (|η| < 1.5) of the MPD experiment have been used. The event plane method uses the correlation of the azimuthal angle φ of each particle with the azimuthal angle Ψ n of the event plane reconstructed from the anisotropic flow itself [2]. The event flow vector for elliptic flow (Q 2 ) and the azimuthal angle of the event plane Ψ 2,TPC can be defined as: where the sum runs over all particles i used in the event plane calculation, and ϕ i and ω i are the laboratory azimuthal angle and the weight for the particle i. The event plane angle Ψ 2,TPC can be used to estimate the magnitude of the elliptic flow v 2 {EP} signal as follows: where R 2 (Ψ 2,TPC ) represents the event plane resolution factor. The cumulants c 2 {k} can be expressed in terms of the moments of the magnitude of the corresponding flow vector Q n ≡ ∑ M i exp(inϕ i ), where M denotes the multiplicity of selected particles in each event [32]. The single-event average two-and four-particle azimuthal correlations can be expressed as follows [32]: For elliptic flow (n = 2), the two-and four-particle cumulants, and the v 2 estimators can be formulated as follows: where the double brackets denote the weighted average of multi-particle correlations over all events.
The non-flow effects may affect the results of the v 2 measurements. They are mainly due to the following particle correlations, and not associated with the reaction plane: Bose-Einstein correlations, resonance decays, and momentum conservation. The multiparticle cumulant removes the contribution of non-flow correlations from lower-order correlations [2] and the v 2 {4} results are expected to be less affected by non-flow effects. In order to suppress non-flow effects in two particle correlation methods: v 2 {2} and v 2 {EP}, one needs to apply the η-gap (∆η > 0.1) between the two sub-events (see [33] for the details).
Elliptic flow fluctuates from event to event and the magnitude of v 2 fluctuations Here, the resulting flow signal, averaged over all events, is denoted as v 2 . In the case of the Q-cumulants (v 2 {2} and v 2 {4}), for a Gaussian model of fluctuations and in the limit of σ v2 v 2 , one can write [2]: This facilitates investigations of the relative fluctuations of v 2 by the lowest ratio of 1, while a weak one leads to v 2 {4}/v 2 {2} ∼ 1. The eccentricity fluctuations make v 2 in the participant plane larger than in the reaction plane v 2 {Ψ 2,TPC } v 2 + 0.5 · σ 2 v2 / v 2 . Figure 8 shows the anticipated performance of the MPD experiment for the v 2 (p T ) measurements of protons and charged pions in Au+Au collisions at √ s NN = 7.7 GeV (upper panels) and √ s NN = 11.5 GeV (lower panels) obtained by four-particle cumulants v 2 {4} (left), two-particle cumulants v 2 {2} (middle) and TPC event plane v 2 {Ψ 2,TPC } (right). The agreement between v 2 (p T derived from a fully reconstructed data analysis based on GEANT4 (open symbols) and UrQMD model data (filled symbols) indicate good performance of the MPD for the detailed differential measurements of v 2 [31,33]. (a,f): four-particle cumulants, two-particle cumulants (b,g) and TPC event plane (c,h). The open symbols correspond to the reconstructed data and closed symbols to the UrQMD model data. for centrality procedure based on the multiplicity of produced particles using the Γ-fit approach are found to be in a good agreement, excluding the 0-10% of central collisions.

Results and Discussion
Here, v 2 (impact) < v 2 (mult) by 4-5%. For the MC Glauber approach, the difference between v 2 (impact) and v 2 (mult) is larger (up to 8-10%) and it has a strong centrality dependence. See panel (b) in Figures 9 and 10 for more information. The larger difference between v 2 (impact) and v 2 (mult) for MC Glauber approach results from the systematic shift between the b and the model data. See Section 3 for more information. The difference between v 2 (impact) and v 2 (mult) does not depend on the flow measurement method and is the same for models UrQMD and AMPT SM. Figures 11-13 show the centrality dependence of v 2 of inclusive charged hadrons from the UrQMD model for Au+Au collisions at √ s NN = 5 GeV (left), 7.7 GeV (center) and 11.5 GeV (right). The results are presented for two particle cumulant v 2 {2} (circles), four particle cumulant

Conclusions
In summary, we have studied the effects of different methods used for the centrality selection on the elliptic flow measurements in Au+Au collisions at √ s NN = 5 GeV, 7.7 and 11.5 GeV within the framework of the cascade version of UrQMD and string melting version of AMPT-SM models. The centralities were defined by the charged-particle multiplicities of produced particles. The elliptic flow v 2 (mult) results from events with centrality classes defined using the MC Glauber and Γ-fit approaches have been compared to v 2 (impact) with centrality classes based on the true impact parameter from the models. The difference between v 2 (impact) and v 2 (mult) is around 1-2% for the Γ-fit approach, except for central collisions where the difference is 4-5%. For the MC Glauber approach, the difference between v 2 (impact) and v 2 (mult) is larger (up to 8-10%) and it has a strong centrality dependence. The v 2 difference increases with decreasing collision energy from √ s NN = 11.5 GeV to 5 GeV for all methods of elliptic flow measurements, used in the present work. The results indicate that the data-driven and model independent Γ-fit approach provides a more accurate way to reconstruct the impact parameter than the model-dependent MC Glauber method. As a result, the centrality selection by the Γ-fit approach will only slightly affect the results of v 2 measurements in 0-10% of central collisions, where the v 2 signal is very small and strongly depends on centrality. Our work may serve as a baseline for the centrality selection of the elliptic flow in the future relativistic heavy-ion collision experiments at NICA energies. In the future, we plan to include the other estimators of centrality based on spectator fragments and extend the study to other flow harmonics.