{'fraction_non_alphanumeric': 0.033693933294997126, 'fraction_numerical': 0.014412018401380103, 'mean_word_length': 4.380874105588861, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generalized Effective String Rope Model for the initial stages of Ultra-Relativistic Heavy Ion Collisions

We present a Generalized Effective String Rope Model (GESRM) for the description of the initial state of relativistic heavy ion collisions. We start from the Effective String Rope Model (ESRM) and take into account fluctuations in the initial state following the Glauber Monte Carlo approach. Results from symmetric Au+Au collisions at different impact parameters and asymmetric A+Au head on collisions are presented at RHIC energies. The produced initial state is used as an initial condition for further hydrodynamical calculations.


I. INTRODUCTION
Relativistic heavy ion collisions allow one to create ultra dense and hot systems, which can undergo a confinement-deconfinement phase transition, leading to the creation of a new state of matter called quark-gluon plasma (QGP) [1][2][3]. In the last two decades with the construction of high energy colliders, such as the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN, as well as with the development of new detectors and new data storing and analyzing methods, the study of heavy ion collisions eventby-event has become possible.
With growing collision energy, the amount of produced particles per event has increased considerably. For example, for central collisions at LHC energies, about 10 4 particles are produced in a single event. This allows sufficiently precise statistical event-by-event analysis of the fluctuations of observables, such as charged particle multiplicity, particle species ratios, transverse momentum, etc. The analysis of these observables would be able to reveal important information about the properties of the system such as, for example, the order of the phase transition, the presence of a critical point, etc. Thus, the study of fluctuations provides us information, which would be unavailable in the study of averages over a large statistical sample of events.
An important observable in relativistic heavy ion collisions is the azimuthal distribution of the emitted particles. There is a principal difference whether the analysis is done on event-by-event basis or averaging over many events. For example, in head-on collisions, i.e. at zero impact parameter, the average overlapping region between the two colliding nuclei possesses azimuthal symmetry while at non-zero impact parameter it has an "almond shape" giving rise to azimuthal asymmetry (see left and center top plots of Fig. 1). On an event-by-event basis, however, the overlapping region fluctuates around the average geometry, leading to azimuthal anisotropy even at zero impact parameter.
The azimuthal anisotropy of emitted particles can be quantified by the Fourier expansion of the particles azimuthal distribution as [4]: where N is the multiplicity of the produced particles of a given type, φ is the azimuthal angle of these particles and ψ n defines the nth-order event plane. For a smooth matter distribution, the second-order event plane ψ 2 coincides with the reaction plane (the plane defined by the beam axis and the impact parameter vector). The coefficients v n are the so-called flow coefficients: v 0 is the radial flow, v 1 is the direct flow, v 2 is the elliptic flow, v 3 is the triangular flow, etc. In Fig. 1  the initial configuration of nucleons of the colliding nuclei [5][6][7][8].
In relativistic heavy ion collisions at non-zero impact parameter, a gradient in the longitudinal component of the flow velocity is present along the impact parameter direction (see right-bottom plot of Fig. 1), the so-called shear flow in z-direction. The largest values of the flow velocity are reached at the extremes of the system, i.e. close to the projectile and target spectators. Such an initial flow velocity configuration has to do with an extremely high angular momentum, which is present in noncentral heavy ion collisions. Due to angular momentum conservation, in the further evolution of the system this may lead to overall rotation of the reaction volume [7] or/and to a large vorticity of the collective flow, which can manifest itself via polarization of emitted particles [9][10][11][12][13][14][15]. This is supported by experiments like STAR collaboration, which has reported observations of global polarization of Λ(Λ) hyperons at non-zero impact parameter in Au+Au collisions [16].
In sec. II we discuss the initial state of ultra-relativistic heavy ion collisions and review the Effective String Rope Model (ESRM) developed to describe such state. In sec. III we present the Generalized Effective String Rope Model (GESRM) and explain how the fluctuations have been implemented into ESRM. In Sec. IV we present the results obtained from the GESRM simulations for symmetric Au+Au collisions at different impact parameters and asymmetric A+Au head on collisions at RHIC energies. We compare the results obtained on event-by-event basis with those obtained averaging over some number of events, N events (N events = 10 − 500.000). Finally, in Sec. VI we summarize and discuss the main results.

II. INITIAL STATE OF ULTRA-RELATIVISTIC HEAVY ION COLLISIONS
The evolution of the system in relativistic heavy ion collisions is commonly divided into three different stages: an initial stage or pre-equilibrium state, an intermediate stage and a final stage or freeze-out.
The initial state describes the first moments of the collision, i.e. from the time when colliding nuclei pass through each other till the time when local equilibration (or at least some pressure isotropization) is established. The equilibration is achieved through the collisions among the constituents of the fireball produced in the initial hard parton collisions. Extreme temperatures and densities generated inside the fireball in ultrarelativistic heavy ion collisions lead to formation of QGP.
An intermediate stage describes the evolution of the fireball from the initial thermalized QGP until the freezeout stage, and nowadays it is usually simulated within a relativistic hydrodynamical model. The outward thermal pressure of the thermalized QGP acts against the inward pressure exerted by the QCD vacuum. Due to the resulting pressure gradients the fireball will undergo 3 dimensional hydrodynamic expansion. As it expands, the temperature decreases and when it drops below a certain value, the so-called critical temperature T c , a deconfinement-confinement phase transition will take place and quarks and gluons will hadronize. The system will further expand and cool down. Finally, the average distance between hadrons will be larger than the strong interaction range, the number of collisions will drastically decrease and finally hadrons will freely move to the detectors. This process is called freeze-out. At this stage a hydrodynamic description of the system is no longer valid and some transport model, such as the Ultra-relativistic Quantum Molecular Dynamics model (UrQMD) [17], should be used to describe the evolution of hadron spectra during the freeze-out process. One of the extreme, but frequently used assumptions, is that freeze out happens on an infinitely narrow hypersurface in the space-time. On its inner part we have a collective matter, described by hydrodynamics, on its outer part we have an ideal gas of different hadron species with momentum distributions generated according to the Cooper-Frye formula [18].
As we can see, relativistic fluid dynamics is commonly used to describe the intermediate stage of the reaction, in which the QGP is assumed to be in a local thermal equilibrium. The final stage of the reaction is simulated either based on freeze out hypersurface or using the hadron cascade model UrQMD. The initial non-equilibrium stage of the collision is, however, the most problematic to be simulated. The main differences among the models are in selecting initial conditions, when the hydrodynamical description becomes valid.
On the other hand we can see that the majority of the most recent models use event-by-event fluctuating initial state. This allows them to analyze simulated events in the same way as the real experiments are performed, as it was discussed in the Introduction, and thus maximize the obtained information. Observing such a trend we decided to update the ESRM for the initial state [20,22], by combining it with the Glauber Monte Carlo approach.

A. The Effective String Rope Model
More than 20 years ago, at the beginning of the RHIC era, a new nucleus-nucleus initial state was proposed, based on the longitudinal Effective String Rope Model [20][21][22] for further 3+1D relativistic fluid dynamical evolution. The big advantage of this initial state, in comparison to others available at that time, was that it reflected correctly not only the energy-momentum, but also angular momentum conservation laws. Consequently, such an initial state for non-central ultra-relativistic heavy ion collisions showed a rather large flow vorticity [9][10][11], and even an effective rotation of the whole fireball has been observed once the ESRM was applied to simulate Pb+Pb collisions at ALICE@LHC [7].
A first step in this model is the creation of a grid in the plane transverse to the beam ([xy]-plane). The collision between two nuclei is then described as a set of independent streak-streak collisions corresponding to the same transverse coordinates {x i , y i }. In the top plot of Fig.  2 such a collision is sketched in the reaction plane (the plane formed by the impact parameter b and the beam direction). The length of each colliding streak is calculated assuming a uniform nuclear matter distribution. Those streaks corresponding to the target, right nucleus, are labeled by l 1 and those corresponding to the projectile, left nucleus, are labeled by l 2 . The space-time evolution of the streaks is governed by the chromo-electric field generated by the color exchange of colliding partons, and by the energy-momentum and baryon charge conservation laws. The bottom plot of Fig. 2 illustrates such an evolution corresponding to two peripheral colliding streaks. At t = 0 both streaks come into contact. The kinetic energy of partons is so high that it is assumed that these can go through the opposite slab of matter without stopping. Only when they have completely passed through each other the chromo-electric flux tube or color string is created (extending up to the target and projectile streak ends), which will slow down and stop the matter. This field is assumed to be uniform as a string, or it is better to say several parallel strings ≡ string rope, because it's string tension, σ, is much higher than that for ordinary hadronic strings. In the ESRM [20][21][22] this uniform field strength was calculated in the following way: where M n is the nucleon mass, ε 0 is the initial energy per nucleon in the Center-of-Mass frame (=Lab frame for RICH), N 1 and N 2 are the baryon charges of the corresponding colliding streaks, and ∆x∆y is the cross section of the streaks. The typical values of dimensionless parameter A are around 0.06-0.08 (for σ measured in GeV /f m). The typical values of σ are 4-10 GeV /f m for ε 0 = 65 GeV, and σ ≈ 6−15 GeV/fm for ε 0 = 100 GeV. We let partons from the colliding streaks evolve according to trajectories derived in ESRM, see [20,22] for more details. They will keep going in the initial direction gradually slowing down up to the some time t = t i,turn , when they will stop, and then they will turn around and go backwards until the two streaks again penetrate through each other, i.e the typical Yo-Yo motion known in the string theory. This scenario of the contraction of the reaction volume is unlikely to occur due to string fragmentations and string-string interactions. In the ESRM it is assumed that at t streak =Max{t i,turn } a uniform streak of length ∆l f , defined by the motion of the outer edges of the colliding streaks, is formed. The uniform energy, e f , and baryon, n f , densities of this final streak as well as its unique rapidity, y f , are calculated from the energy, momentum and baryon charge conservation laws, and the corresponding Equation of State (EoS). We assume for the QGP a Stefan-Boltzmann gas EoS, i.e. p = e/3 − 4B/3, where p is the gas pressure density, e the energy density and B the Bag constant (more details will be shown in section IV C).
Once this final streak is formed it starts to expand into the vacuum with velocity of light according to the analytical solution [22,23].
The bottom plot of Fig. 2 shows the trajectory of leading partons from t = 0 up to t = t streak , blue solid lines, and the further expansion, red dashed lines. The analytical solution for the expansion of the homogeneous final streak into the vacuum [22,23] allows us to know the energy and baryon density distributions and the local rapidity along the streak. This is valid until the expansion waves coming from the edges with the velocity of sound did not cross in the middle of the streak.
Thus, considering the overall reaction volume, we can present the initial state, i.e. energy density, baryon density and flow distributions of the fireball at some time moment t = t f in . This time, t f in is well chosen if at this moment the final streaks have already been formed and started their expansion for most of the transverse coordinates, and, at the same time, none of them has reached already the moment when the analytical expanding solutions are not applicable [22]. Please note that by assumption we neglect the transverse expansion in the model until t = t f in ; thus, at this moment all the cells have only longitudinal velocity.
In Figs. 3 and 4 we show the energy density and velocity distributions in the reaction plane obtained from the ESRM with final expanding streaks for Au+Au simulations at initial energy ε 0 =100 GeV per nucleon (i.e. √ S N N = 200 GeV) for impact parameter b 0 = 0.5. In general case in our notation the module of the impact parameter b is given by From the energy density distribution one can see that, at non-zero impact parameter, the system forms a type of tilted disk, and, thus, the direction of fastest expansion, generated by the strongest pressure gradient, will deviate from both the beam direction and the transverse flow direction, giving rise to a new flow component called third flow component or antiflow [24].
The velocity distribution shows that while in the central zone of the collision partons move rather slowly, those at the ends of the streaks move with a much higher velocity. Obviously the initial flow distribution, shown in Fig. 4, will generate a high vorticity; this quantity will be discussed in more detail in section IV D.

III. GENERALIZED EFFECTIVE STRING ROPE MODEL
In this section we present the Generalized Effective String Rope Model, in which the fluctuations in the initial state of relativistic heavy ion collisions are taken into account following the Glauber Monte Carlo approach [25].  [22,23]. Finally, at t = t f in = 5 fm, our final streak has an energy density profile shown at the top of the plot.

A. The Glauber Monte Carlo approach
In the Glauber Monte Carlo approach fluctuations are introduced randomly distributing positions of the nucleons, which allows us to obtain different configurations of those for each colliding nucleus in each collision. Therefore, the number of participant nucleons will now fluctuate event-by-event leading to fluctuations in the different physics quantities such as baryon charge, energy, the total momentum and consequently the central rapidity. numbers {x, y, z, δ}, where {x, y, z} are Cartesian coordinates from 0 to R, and δ is a random number, which takes the values in the range 0 ≤ δ ≤ 1. If the following condition is satisfied then we take the numbers {x, y, z} as the coordinates of the center of a given nucleon, otherwise, we generate a new set of random numbers until this condition is satisfied. In this way we ensure that, in the average, our nucleons are distributed according to Woods-Saxon (WS) nuclear matter density distribution. The skin depth a and the nuclear radius parameter R have been taken from Ref. [25] corresponding to a 197 Au nucleus, namely a Au = 0.535 fm and R Au = 6.38 fm. Please note that using such a WS baryon distribution, while in the ESRM colliding Au nuclei have been considered as a homogeneous spheres with the radius of 7 fm, we change a bit the geometry of the collision, and now impact parameter b 0 = 0.5, i.e. b = 0.5(R Au + R Au ) = R Au = 6.38 fm, doesn't mean the same as in original ESRM calculations [20,22].
The second step is to consider each nucleon as a sphere of radius R N = 0.842 fm [26], in which the nucleon baryon charge is homogeously distributed. In Fig. 5 is shown an example of two given configurations of nucleons corresponding to two 197 Au nuclei just before a collision at impact parameter b 0 = 0.5. These distributions are shown in the reaction plane, top plot, and in the transverse plane, bottom plot. As can be seen, the "almond-shape"of the overlapping region between both nuclei is not observed on event-by-event basis. Instead, this will fluctuate around the average geometry producing any order geometric deformations. The almond-shape of the overlapping region can be restored averaging over many events.

Nucleonic matter discretization
Once the nucleons have been randomly distributed, we generate a two-dimensional grid in the transverse plane, in which the cell-size is taken to be: ∆x∆y = (R Au /10)(R Au /10), and we calculate the baryon charge in each cell. However, there is no easy analytical way to perform this calculation. So instead we will follow a Monte Carlo approach and consider each nucleon as a set of fictitious particles randomly distributed within a homogeneous sphere of R N radius. This method allows us to obtain the baryon charge in each transverse cell in an easy and fast way.
Each of these fictitious particles are associated with a certain baryon charge whose value depends on the total number of these. In our case, we have used N =1000 fictitious particles so the baryon charge associated to it is 10 −3 . The total baryon charge in each transverse cell will be then the sum of all fictitious particles within it, and the sum of the nucleons over all transverse cells will give us the total baryon charge of the nucleus, that in the case of 197 Au is 197.
For a given nucleon, we generate a set of three random numbers {x, y, z} within the intervals; where R N is the nucleon radius and {x i , y i , z i } are the coordinates of the center of the ith nucleon. If the set of random numbers {x, y, z} satisfied the relation we take these numbers as the coordinates of a given fictitious particle, otherwise, we generate a new set of random numbers until this condition is satisfied. In Fig. 6 an example of such a distribution for 1000 fictitious particles is shown for a given nucleon. The top plot illustrates the fictitious particle distribution in 3 dimensions, while in the bottom plot the same distribution projected in the transverse plane is represented. The size of each cell of the grid has been chosen to be equal to the one used in our calculations, i.e. ∆x = ∆y = R Au /10 = 0.638 fm.

B. Implementation of fluctuations in the ESRM
In the ESRM the starting point was related with the lengths, l 1 and l 2 , of each colliding streak. Now we have the baryon charge in each transverse cell and its distribution along the z-direction. In the top plot of Fig. 7 an example of the baryon charge distribution in the reaction plane for a nucleus at rest is given. As we can see, due to the random position of nucleons, there are some regions along the z-axis with a high baryon charge densities, and there are others with zero baryon charge. However, due to the strong Lorentz contraction and further collision, the initial z-distribution of the baryons/quarks is not so important. What is important is the number of baryons/quarks on each side, which will collide with each other. Thus, we can assume that the baryon charge is homogeneously distributed in a volume ∆x ∆y l, where ∆x ∆y is the transverse area of each cell and l is the  length of the streak, which is given by where ρ 0 is the normal nuclear matter density and N α i,j is the baryon charge in the transverse cell (i, j) corresponding to the αth nucleus. In the bottom plot of Fig. 7 the length of each streak, obtained from the total baryon charge distribution of the top plot, is illustrated. Due to fluctuations, the largest streaks are not necessarily located in the middle region of the nucleus but these can be now in any place, even at the top and bottom extremes, leading to fluctuations around the average geometry. Thus, the spherical symmetry assumed in a 197 Au nucleus at rest is no longer true on event-by-event basis. However, averaging over many events the geometry can be restored. Describing a collision of two nuclei which are moving with relativistic energies it is necessary to take into account the Lorentz contraction, and, thus, Eq. (6) is generalized to Fluctuations on the number of colliding baryons, which is equivalent, according to Eq. (7), to fluctuations in the lengths of colliding streaks, could lead to formation of holes and cells with low baryon charge in the central zone of the overlapping region between both colliding nuclei. The string tension, σ, of the chromo-electric string formed from the collision of two of such streaks will be very small. On the other hand, to maintain the average, in almost every event there appear some streaks (one or more) which generate a rather high string tension, and consequently will form final streak in a short time and, as it was discussed above in section II A, these will limit the maximal time when we should stop our model and fix the produced initial state.
If the old definition of string tension is used, see Eq.
(2), then, first of all, we have to reduce the final time, t f in , from 5 fm, to 4.5 fm, for A = 0.65, or even less for higher A. At the same time we observe that there was not enough time to start the final streak expansion for most of the transverse cells.
To ensure that most of the final streaks have been formed and have started their expansion a modification of σ is necessary. In the present work we have used the following parametrization for σ: where A is now a dimensional parameter measured in GeV; in the calculations we used A = 0.05/ √ ∆x∆y = 0.0784 GeV. This new expression for σ presents a lower dependence on the product N 1 × N 2 , and, thus, the difference between lower and higher string tension will be much less. In consequence, the final streak formation times for different transverse coordinates will be much more homogeneous. In Fig. 8 both definitions of σ as a function of the product of the baryon charges inside the colliding streaks (N 1 and N 2 ) are presented. The expression used in the ESRM has been indicated by σ old while that used in the present work is indicated by σ new . In both cases, the same values have been used to the parameters present in these expressions: A = 0.0784 (in GeV for the new definition of σ (see Eq. 8)), ε 0 = 100 GeV and ∆x∆y = 0.407 fm 2 .

IV. RESULTS
In this section we present the results obtained from GESRM. We compare those obtained on event-by-event basis with those averaging over N events, where N will vary from 10 to 500.000 events. The latter should be qualitatively comparable with the ESRM results, but we note that even if we average over very large number of events we will not reproduce ESRM case, and this has to do not only with the new definition of σ (Eq. (8)) and WS distribution of nucleons in Au. Actually these modifications would only generate a small quantitative, but not qualitative difference. On the other hand, initial state fluctuations lead to principle differences in the initial state not only for a single (random) event, but also if we perform an averaging over many events. As an illustration of this statement let us consider the final profile of the energy density at t = t f in for some given transverse coordinates. The result of ESRM is well known and shown by a blue dashed line in Fig. 9. In the generalized model the situation will change: the number of colliding nucleons (from both sides) will now fluctuate; this will lead to different field tensions σ and correspondingly different final streak length in each event, as well as to a different rapidity of the final streaks y f (which is defined by the momentum conservation). If we average over many such events then, according to the Central Limiting Theorem, we can expect to see some Guassian-like shape. This is exactly what we see in Fig. 9, where green solid line corresponds to an averaging over N = 10000 events.
Generalizing this discussion for the whole reaction volume is shown in Fig. 10, where the energy density distributions in the reaction plane, obtained from the ESRM, top plot, and from the GESRM averaging over N = 10000 events, bottom plot, are shown for ε 0 = 100 GeV and b 0 = 0.5. In both cases, we used the same values of parameters A = 0.0784 (in GeV for the new definition of σ (see Eq. 8)). The averaged GESRM initial state is a bit wider and less peaked (=more smoothed) in the middle.

A. Reaction volume and number of participant nucleons
We consider the reaction volume as the region of space occupied by our expanding system at t f in = 5 fm. Unlike the ESRM, where the reaction volume will be fixed for a given impact parameter, in the GESRM, due to the random distribution of nucleons, this will fluctuate eventby-event. The bottom plot of Fig. 11 is representing the initial state volume (V IS ) as a function of the impact parameter for N = 1, 10, 100 and 1000 events. For a single event, most of the non-empty cells are located in the overlapping region and only a few are outside (see Fig. 5). However, when we average over many events, the contribution of the outside region becomes more and more significant giving rise to a considerable increment in the reaction volume. Please note that this is just a pure effect of the geometric fluctuations of the reaction volume. This doesn't mean that the number of participants grows, since those peripheral cells have very small densities.
This behavior of the number of participant nucleons (N part ) is clearly seen in the top plot of Fig. 11. To obtain N part we have calculated the baryon charge density n i inside each non-empty cell and we have multiplied it by the volume of a given cell in Lab frame: where V cell = ∆x∆y∆z, N cell is the total number of nonempty cells and γ i is the Lorentz factor of the cell i. We can see that for a single event the fluctuations of N part increase when we increase impact parameter, and reaches its maximum around b 0 = 0.5. On the other hand, with increasing number of events the number of participant nucleons rapidly converge to the average value, which obviously depends on the given impact parameter. Please note that in the generalized model for central collisions (b = 0) the average number of the participants is not N part = 197 + 197 = 394, but, due to initial state fluctuations, it is a bit smaller: In order to check the stability of our calculation we performed simulations with different cell sizes, namely ∆x = ∆y = R Au /6; R Au /10 and R Au /18. The corresponding results are shown in Fig. 12. Both the number of participants, top plot, and the reaction volume, bottom plot, grows a bit for a bigger cell size, what can be expected. We also note that this difference is reduced for the most peripheral collisions.

B. Baryon charge and energy density distributions
The baryon charge and energy density distributions at t = t f in are the main output of the GESRM, which should be used as an input for further hydrodynamical evolution.
In Figs. 13 and 14 we represent the baryon charge and the energy density distributions for a single event with fluctuations, top plots, and averaging over N = 10000 events, bottom plots. These are shown in the reaction plane, left plots, and in the transverse plane, right plots. Similarly to the results obtained from the ESRM, when we average over many events the baryon charge and energy density distributions show a type of tilted disk in the reaction plane. This is not clearly seen on eventby-event basis because of fluctuations. With respect to the transverse plane, we see that for N=10000 events, the overlapping region is "almond shape"giving rise to a strong elliptic flow. However, on event-by-event basis, the geometry of the collision will fluctuate around the averaged geometry, giving rise to geometric deformations: elliptical ε 2 , triangular ε 3 , quadrupole ε 4 and other harmonics, which will generate momentum anisotropies quantified by flow harmonics v n [27] as discussed in Eq. (1).
In Figs. 15 and 16 are shown the longitudinal linear momentum and the rapidity distributions, respectively. As we can see, when we average over N = 10000 events the fluid cells in the central zone of the collision are approximately at rest because the colliding streaks are symmetric in this region. This is not observed on event by event basis because of fluctuations of the streak lengths.
There is a remarkable difference between the distributions obtained averaging over many events, bottom plots of Fig. 16 lision, where the streak ends move in opposite direction as well.
From the Central Limiting Theorem it is expected that fluctuations in the different physics quantities leads to a Gaussian-shape. This is illustrated in Fig. 17, where the distributions of the results from different single events for the total longitudinal linear momentum of the collisions (top plots), its total energy (middle plots), and center of mass rapidity of the system (bottom plots) are shown at impact parameters b = 0 fm (left plots), and b = (R Au + R Au )/2 fm (right plots). At each histogram of Fig. 17 we present 10.000 results. For the b = 0 case, as it is expected for symmetric collisions, the Gaussians corresponding to the total linear momentum and to the center-of-mass rapidity are centered around 0. In the case of the energy, the Gaussian is centered around the average energy of the system corresponding to zero impact parameter, and it can be calculated following the result of Fig. 11: It is interesting to compare the widths of the corresponding Gaussian distributions; these clearly grow with impact parameter. This is not a surprise since we have seen already that the relative fluctuations are stronger for smaller systems. To be more quantitative in Fig.  18 we present the Gaussian width corresponding to the center-of-mass rapidity (δy) as a function of the impact parameter.
The results of Fig. 18 show that in Ref. [7] the authors have overestimated the center-of-mass rapidity fluctuations, taking δy = 1 and 2, even for very peripheral events.

C. Chemical potential, Temperature and Entropy
One can obtain temperature and baryon chemical potential from the energy and baryon charge densities, calculated in our model, applying the corresponding EoS. EoS given by the Stefan-Boltzmann expressions: where T and µ q are the quark temperature and chemical potential (µ b = 3µ q ), and n SB = n b is the baryon charge density in the quark phase. From Eq. (12) it is possible to write the temperature as a function of n SB and µ b , and replace it in Eq. (10) in order to obtain µ b . For N c = N f = 3 we find the following expression: Please note that in the MIT Bag Model EoS, used for our calculations, the energy density and pressure of the parton gas differ from e SB and P SB by the Bag constant B: In our calculations we have used B = 0.330 GeV/fm 3 . If the pressure is high enough then e SB = e − B and we can proceed solving Eq. (13). However, if P SB = e SB /3 < B we might run into the problem of negative pressure. Following the prescription given in Ref. [28], we assign in such a situation P = 0 and e SB = 3B e 4B 4/3 .
Once we have calculated the baryon chemical potential and the temperature, we can easily obtain the entropy density, s, from the following thermodynamical relation: In Figs. 19, 20 and 21 the baryon chemical potential, µ b , the temperature, T , and the entropy density, s, distributions are represented respectively. As we can see, the temperatures reached in the middle region of the collision are very high -up to 450-500 MeV. From the entropy density one can calculate the total entropy per baryon charge, S/N : In Fig. 22 we present the entropy per baryon charge in our initial state, simulated for Au+Au collisions, as a function of the initial energy per nucleon at impact parameter b 0 = 0.5, top plot, and as a function of the impact parameter at √ S N N = 200 GeV, bottom plot. In the top plot we can see that, first of all, S/N increases linearly with ε 0 . This linear rise, observed for a given impact parameter, i.e., given average < N part >, is directly related to the rise of entropy with average temperature of the initial state in more energetic collisions. Secondly, we note that S/N saturates with the number of events, although not very rapidly -the saturation is reached around N events = 100, but the saturated value is much higher than S/N in a single event! For a single event we can also notice that fluctuations are weakly dependent on ε 0 . In Ref. [28] the value of S/N 225 for the initial state calculated in ESRM model for Au+Au collision at ε 0 = 65 GeV per nucleon for central collision, b 0 = 0. This value is compatible with our new results in the generalized model for a single event. However, for the initial state averaged over many events the S/N is much bigger. Thus, although the averaged initial state looks like smoothed over corresponding initial state obtained in ESRM, Fig. 10, it will lead in a subsequent hydrodynamic evolution to rather different results due to much higher initial entropy.
In the top plot of Fig. 11 we have seen that for a given impact parameter, N part saturates very quickly with number of the events. Furthermore, the bottom plot of Fig. 11 shows how the reaction volume increases with the number of events for given impact parameter. We think that this is the effect, that explains the rise of saturated value of S/N with respect to that in a single event. Since N part saturates rapidly the observed fluctuations in S/N are related to fluctuations of the entropy of initial state. To study these fluctuations in more details we present these for different system sizes in the bottom plot of Fig. 22.
In the bottom plot of Fig. 22 one can distinguish two regions: the first region corresponds to 0 ≤ b 0 < 0.6, where the S/N is practically independent on b 0 , and the second one with b 0 ≥ 0.6, where a strong dependence on b 0 and on the number of events is observed. The more peripheral the collision is, the greater is the entropy generated by our model per baryon charge. Please note that also the fluctuations of S/N for a single event strongly increase with impact parameter. To understand these results we have to recall the basic thermodynamic relation: To estimate the average variation of S/N we first calculate the corresponding differential and using Eq. (17) we finally obtain the following relation From Fig. 11 we know that the fluctuations of relative volume and of N part increase for smaller systems, i.e. for higher impact parameters, and thus, following Eq. (19), the S/N fluctuations in a single event will increase with b 0 , as it is seen in the bottom plot of Fig. 22. However, why the S/N value itself grows with impact parameter for peripheral collisions? As impact parameter grows, ∆b 0 > 0, the number of participants and the reaction volume decrease, ∆N < 0 and ∆V < 0, and the average change of S/N can be estimated from Eq. (19): (20) Since P T > 0 the second term always generates a decrease of S/N . On the other hand the first term can be both positive or negative, and for our simulations it is positive. As we can see in the top plot of Fig. 23 the first coefficient 100 GeV T − S N < 0 for all the values of the impact parameter, and, thus, the first term is always positive. However, in the first region (0 ≤ b 0 < 0.6) this increase of S/N is compensated by the second negative term, see the bottom plot of Fig. 23. As we go to more and more peripheral collisions the average temperature slowly drops down while S/N slowly grows maintaining the first term practically constant, while the value of the second term decreases, since P T ∼ T 3 . Thus, at some moment (b 0 0.6) the S/N starts to grow and then it grows faster and faster, just as we observe in Fig. 22, since according to Eq. (20) this is an auto-reinforced process. The bigger S/N is, the bigger is the first positive term and, thus, the bigger is the increase of S/N . We have also checked whether the S/N production is sensitive to the cell size. Figure 24 shows that it has a rather little sensitivity. The most central events (b 0 ≤ 0.5) are practically independent to the cell size, and for more peripheral collisions we observe a small increase of S/N with decrease of cell size. Taking into account that from Fig. 12 we know that N part for very peripheral collisions (b 0 ≥ 0.9) is insensitive to the cell size, Figure  24 tells us that the entropy production in the peripheral collisions grows a bit if one considers smaller cells.
We can also see in Fig. 24 that the important effect, seen earlier in the bottom plot of Fig. 22, that entropy per nucleon production for the initial state averaged over many events strongly increases with respect to that for a single event, does not depend on the cell size.

D. Vorticity
The vorticity field is a quantitative measure of the local circulation of a fluid. It is calculated in each point of the fluid and can be shown via the so-called vortex lines. A vortex line is a line whose tangent is everywhere parallel to the local vorticity vector. Vortex lines cannot cross each other, which implies that vortex lines cannot emerge or terminate anywhere of the fluid, but must keep on going until they reach the boundaries of the flow.
In classical physics for incompressible, perfect fluids vorticity exhibits an impressive conservation law: the conservation of circulation [30]. The relativistic fluid dynamical calculations indicate [29] that typical flow patterns and instabilities may occur here also. Thus, their studies can provide insight into the properties of the QGP.
In peripheral heavy-ion collisions the initial angular momentum of the system generates a strong shear and vorticity in the flow [9], which may lead to rotation [7] and even Kelvin Helmholtz instability [29] in the reaction plane for low-viscosity QGP. At later stages of the reaction particles produced in the vortical matter are expected to be polarized.

Classical vorticity
Mathematically, the classical vorticity ω is defined as the curl of the flow velocity v: Thus, for example, in the reaction plane the vorticity will be where the x,y,z components of the 3-velocity v are denoted respectively. In this definition of vorticity we have already included the factor 1 2 for the symmetrization to have the same magnitude of vorticity as for symmetrized volume divergence or expansion rate [31].
To calculate the y-component of vorticity, ω y , in each cell of the reaction plane we proceed as follows: First we label each cell by the indexes i,k corresponding to x,z axes, respectively. For a given layer y, we have a contribution from the side points 1,2,3,4 and the corner neighboring points 5,6,7,8 (see Fig. 25).  Since in our model there is only expansion in the zdirection, the first term on the right side of Eq. (22) vanishes. From this and the definition of partial derivative we can rewrite Eq. (22) in the form We now take all possible differences between neighbors cells in the x-direction, i.e., we take the differences: 1-0,5-2,2-6,0-3,4-7,8-4. So the vorticity for a given i,k cell will be given by This is not true for surface cells since we must delete every term where there is an empty cell, since to replace the empty cell by zero velocity would lead to a large derivative. Obviously the process is the same for all ylayers. In Fig. 26 the classical vorticity distributions obtained for a single event, top plots, and averaged over N =50000 events, bottom plots, are shown. Please note that, from the angular momentum conservation, the overall vorticity ω class y should be negative, what is clearly seen for the averaged initial state. We can note that for N =50000 events vorticity takes positive values only at the extremes of the reaction volume, as it has been seen in the other smooth initial state models, for example [9,32], while for a single event positive values can be distributed throughout all the reaction volume due to fluctuations. If we look at the average vorticity in the transverse plane, we can see two well separated regions: a central one, where the vorticity takes small and negative values, and an outer one, where strong fluctuations are observed. However, the cells corresponding to this outer region contain small amount of matter and energy (see Figs. 13 and 14). So it is expected that if one studies the energy or matter weighted vorticity these strong fluctuations will vanish.
In Fig. 27 the results obtained for the y-component of the classical vorticity summed over all y-layers are presented, i.e.
The top plot shows Ω class y for a particular single event with fluctuations, while the bottom plot presents Ω y averaged over N =50000 events. Again, in the case of Au+Au, √S NN = 200 GeV, b=(R Au +R Au )/2, t=5 fm

Relativistic vorticity
Unlike in classical hydrodynamics, where vorticity is defined as the curl of the velocity field v, several vorticities can be defined in relativistic hydrodynamics [33]: the kinematical vorticity, the kinematical transverse vorticity, the T-vorticity, the thermal vorticity,... In this work we follow the definition used in Ref. [31], known as the transverse kinematical vorticity for the relativistic case. The corresponding vorticity tensor is defined as where for any 4-vector q µ the quantity ∇ α q µ ≡ ∆ β α ∂ β q µ = ∆ β α q µ ,β and ∆ µν ≡ g µν − u µ u ν . This leads to where ∂ τ u µ ≡u µ = u α ∂ α u µ is the proper time derivative of u µ .
From the rapidity distributions represented in Fig. 16 we can see that the rapidities in the central region of the collision are rather small. This indicates a strong stopping of partons in the middle region of the reaction volume and, therefore, it is expected that in this region the acceleration of the fluid elements are negligible compare to the rotation: |∂ τ u µ | << |∂ x u z |. This is also true at the extremes of the system where partons move with almost no loss of rapidity. This approximation allows us to simplify the relativistic vorticity to the following expression Thus, for example, for vorticity development in the reaction plane we have Since for each streak the cells moving at velocities close to the speed of light are those located at the very ends of each streak, it is expected that these cells will have a much larger vorticity in relativistic case than in the classical one. In Fig. 28 the y-component of the relativistic vorticity (ω rel y ≡ ω z x ) is shown. If we compare classical and relativistic vorticity distributions in the reaction plane, we will see that the latter at the edges of the reaction volume (at large |z|) is one order of magnitude bigger than the former because of the gamma factor present in the relativistic equations. On the other hand, in the center of the reaction volume both vorticities have similar values, since in this region the gamma factor is close to unity.
Of course, the cells showing the biggest relativistic vorticity contain a very little amount of matter, and the situation might drastically change if one studies an energy or matter weighted relativistic vorticity, but we are not going to enter here in such a discussion.
In Fig. 29 we show the relativistic vorticity summed over all y-layers (Ω rel y = i ω rel yi ) in the reaction plane for a single event, top plot, and averaged over N =50000 events, bottom plot.
In Fig. 30 we present the classical, top plot, and relativistic, bottom plot, total vorticities, i.e. sum of the corresponding vorticities ω class/rel y over all non-empty cells of the grid: Au+Au, √S NN = 200 GeV, b=(R Au +R Au )/2, t=5 fm as a function of the impact parameter for 1, 1000 and 50000 events. As we can see, the results show that the total vorticity is strongly dependent on the impact parameter, similar to the angular momentum (see for example [34]), presenting a parabolic-like shape with a shallow minimum located around b 0 =0.75, for the classical total vorticity, and b 0 =1.0, for the relativistic total vortic-ity (please note that both, relativistic and classical, total vorticities show approximately the same shape, although they differ in their magnitudes by factor 2).
One interesting thing that can be also noticed from Fig. 30 is that the total vorticities show very strong fluctuations even when we average over N=50000 events. In fact, there is a rather small difference between the spread of the points obtained for N=1000 events and those obtained for N=50000 events, what suggests a very slow "saturation" of the total vorticity with the number of events.
This is actually an effect that can be expected. Making an initial state averaged over some number of event, what we actually average is the baryon and energy content of each cell, assigning to it < N 0 >, < N z >, < T 00 > and < T 0z >. These quantities, directly related to the corresponding conservation laws, relatively quickly saturate to their medium values, determined by the impact parameter, so as the quantities directly related to them, like, for example, flow velocity < v z >=< N z > / < N 0 >. However the vorticity has to do not with the flow velocity, but with its derivative, which naturally will need much more statistics to demonstrate the saturation effect. To understand better such large fluctuations of vorticity, we show in Fig. 31 distributions of the classical vorticity single event results for Au+Au collisions at √ S N N =200 GeV, t f in = 5 fm for b = 0 (top plot) and b = (R Au + R Au )/2 (bottom plot). The total number of results presented at each histogram is 500.000. As we see, for such a large number of events we obtained a Gaussian-shaped curve centered at ω class = 0 fm −1 , for b 0 = 0, and ω class ∼ -65000 fm −1 , for b 0 = 1/2, as it should be according to Central Limiting Theorem. However we can also note that the standard deviation in the case b 0 = 0 is surprisingly high -σ = 38036±68 fm −1 . For b 0 = 1/2 it is approximately two times smaller -σ = 18123±38 fm −1 , but still the distribution is rather wide, if we compare for example with those of Fig. 17. Out of this we can guess that in order to get some experimental results related to the vorticity with good accuracy a rather high statistics may be needed. Looking at Fig. 30 we can also note that the total vorticity for an averaged state may be different from the vorticity in a single event. This can be seen despite the strong fluctuations. This is the same effect as the one seen already for the entropy per baryon production, and, in our opinion, it has the same reason. From Figs. 28 and 26 we have learned that the most peripheral cells have the highest vorticity, and not always negative! In the IS averaged over many events the reaction volume grows, see Fig. 12 (bottom plot), due to an increment of peripheral cells with a little amount of matter, but with rather high vorticity, what leads to the observed results. On the other hand since these cells contain a little amount of matter then their contribution, for example, to Λ hyperon polarization results will be very little, since the probability that such a heavy baryon will be produced in one of these very peripheral cells is very low.
Finally, we decided to test how sensitive are the obtained vorticity values to the cell size of our model. In Fig. 32 we show the classical, top plot, and relativistic, bottom plot, vorticities summed over all non-empty cells of the grid and normalized to total number of cells (which will obviously vary with the cell size), as a function of the impact parameter for different cell sizes. Despite the fluctuation we can clearly see that the generated vorticity, averaged over 50000 events, is not sensitive to the size of numerical grid (not more than any other derivative).

E. Initial state for asymmetric A+Au collisions
In order to see more clearly which features of the initial state are sensitive to the geometry of collision and which depend only on the number of participants we compare peripheral Au+Au with asymmetric central A+Au collisions for different A from 2 to 180. The idea is to compare initial states generated in these geometrically rather different situation but with the same average number of participants. The corresponding initial state calculations are performed for √ S N N = 200 GeV and t f in =5 fm; the average number of participants is compared for N events = 100.
To parameterize the nucleon density distribution of the nuclei with A nucleons we have used: -for A≥100, a WS distribution with parameter a W S = 0.535 fm and radii calculated from equation r = 1.1A 1/3 fm.
-for 20≤A≤100, a WS distribution with parameter a W S = 0.535 fm and radii calculated from equation r = 1.2A 1/3 fm.
-for 6≤A≤20, an Harmonic Oscillator (HO) distribution where the parameters α and a are presented in Tab. I. -for A=2, 3, 4 we have used an homogeneous distribution with radii R 2 = 2.095 fm, R 3 = 1.976 fm, R 4 = 1.696 fm correspondingly, taken from Ref. [35]. In order to be able to compare the results for peripheral Au+Au with asymmetric central A+Au collisions, we, first of all, perform a polynomial fit to the average number of participants for Au+Au collisions, presented in Tab. II, as a function of the corresponding impact parameter. The resulting fit is shown in top plot of Now, when we know that the average number of participants in the corresponding peripheral Au+Au and central A+Au collisions is the same, we can look at the reaction volume and the entropy per nucleon of these initial states. The corresponding quantities are shown in Fig.  34. As we can see, for both these observables the increase of the impact parameter in symmetric reaction is equivalent to the reduction of the projectile size in asymmetric reaction for central and semicentral collisions, i. e. for b 0 ≤ 0.7 or in other words for A ≥ 20. And this equivalence is true not only for the averaged over 100 events initial state, but also for single events.
For a very peripheral (lowest A) collisions we do observe some difference. For the reaction volume this difference in not very significant. The fact that for the three smallest A the reaction volume does not grow when we average over 100 events has to do with the fact that for these nuclei their A nucleons are distributed homogeneously within a sphere of a rather small radius. Thus, the fluctuations of initial positions of their nucleons can not increase the maximal volume of the fireball, as this happens for the distributions which do not have a well defined limits.
On the other hand, the entropy per nucleon production in a very peripheral symmetric collisions is principally different from that in central asymmetric collision with the average number of participants (since N is the same, all the observed difference in S/N has to be attributed to the entropy production). In the peripheral Au+Au collisions S/N shows a very strong rise, which we have seen and discussed in section IV C (see Fig. 22), while in asymmetric central collisions of small nuclei with Au the S/N remains more or less constant.
We can also see in Fig. 34 again an important effect  that the S/N production for the initial state averaged over many events strongly increases with respect to that for a single event, similarly as it appear in the bottom plot of Fig. 22. However, if we compare the classical and relativistic vorticities (summed over all the cells), which are shown in Fig. 35, we will observe that the results obtained in asymmetric A+Au head on collisions are completely different from those obtained in symmetric Au+Au collisions at different impact parameters. While in the latter the vorticity decreases with impact parameter until it reaches some minimum value, in the central A+Au collisions it oscillates around zero, as it should be due to symmetry reasons. Again, as in the previous section, we observe a qualitatively correct (expected) behavior of the vorticity accompanied by very strong fluctuations.

V. SUMMARY AND DISCUSSIONS
We have presented a generalized ESRM, which takes into account fluctuations in the initial state of relativistic heavy ion collisions following the Glauber Monte Carlo approach. The random distribution of nucleons leads to fluctuations in the geometry of the collision. In particular, the overlapping region of the two colliding nuclei shown in Fig. 5 for a single event deviates from that obtained for homogeneous spherical nuclei. The initial state fluctuations lead to an enhancement of the different components of the azimuthal distribution of emitted particles; the odd components, such as direct flow (v 1 ) and triangular flow (v 3 ), will be particularly affected.
By averaging our initial state over large number of the events we observed that the average values of the conserved quantities, like N part , rapidly converge to some value, determined by the initial energy and impact parameter of the simulated reaction; and these values can be compared with the corresponding value in the initial state model without randomization. Qualitatively such a comparison is quite fair, although they never will be identical. For example, even for the N part , which converges very rapidly, there is some discrepancy -< N part (b = 0) > will never be equal a A 1 + A 2 , as it was shown in section IV A. The energy, baryon charge and flow velocity distributions generated in GESRM, averaged over many events, also look qualitatively similar to those of ESRM, although the overall volume of the fireball is much bigger -such an averaging adds many cells with a very little amount of matter.
On the other hand we have shown in section IV C that there is a principle difference in using single event initial state and the one averaged over many events with randomization. This difference mainly has to do with an entropy per nucleon, which strongly grows in an averaged initial state. Similar trend has also been observed for asymmetric A+Au collisions. We would like to stress that the S/N in a single event in generalized ESRM with randomization and in the original ESRM, considering colliding nuclei as homogeneous spheres, are fairy comparable. Thus, GESRM initial state, averaged over many events, only looks like smoothed over corresponding initial state obtained in ESRM (Fig. 10) -due to much higher initial entropy it will lead to differences in hadron production.
Generalized ESRM, as well as ESRM, generates a strong sheer and vorticity in the initial state. We have seen in section IV D that in our model the total vorticity shows very strong fluctuations. These strong fluctuations require averaging over N 50.000 events in order to see any sign of saturation. However, the presented results