{'fraction_non_alphanumeric': 0.03962365180142278, 'fraction_numerical': 0.01455289528034881, 'mean_word_length': 4.144923258559622, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large Baryon Densities Achievable in High Energy Heavy Ion Collisions Outside the Central Rapidity Region

Nuclei are nearly transparent to each other when they collide at high energy, but the collisions do produce high energy density matter in the central rapidity region where most experimental measurements are made. What happens to the receding nuclear fireballs? We calculate the energy loss of the nuclei using the color glass condensate model. We then use a simple space-time picture of the collision to calculate the baryon and energy densities of the receding fireballs. For central collisions of large nuclei at the BNL Relativistic Heavy Ion Collider and the CERN Large Hadron Collider we find baryon densities more than ten times that of normal nuclear matter. These results provide initial conditions for subsequent hydrodynamic evolution and could test the equation of state at very high baryon densities.


I. INTRODUCTION
Since the 1980's the high energy heavy ion community has focussed on the central rapidity region of nucleus-nucleus collisions [1]. That followed the influential work of Bjorken in 1983 [2]. The reasons are: (i) the energy density is expected to be higher there; (ii) the matter is nearly baryon-free, making it relevant to the type of matter that existed in the early universe and more amenable to comparisons with lattice gauge calculations; (iii) detectors in a collider can more readily measure particle production and correlations within one unit of rapidity around the centerof-momentum of the colliding nuclei. The earlier work of Anishetty, Koehler and McLerran in 1980 [3], which found that nuclei were significantly compressed and excited when they collide at extreme relativistic energies, was pursued only sporadically. In Ref. [4] we picked up the line of work of Ref. [3] in light of the theoretical and experimental developments since then. We found that in central collisions of gold nuclei at √ s N N = 200 GeV, the top RHIC energy, the central core lost about ∆y = 2.9 units of rapidity and was compressed by a factor of e ∆y [3][4][5][6], a very large compression indeed. In this paper we expand our work on the baryonic fireballs which emerge beyond one unit of rapidity. It should be noted immediately that our calculations are not particularly relevant to the lower beam energy scan at RHIC.
In our work we use the semi-analytical results of Refs. [7,8], which were based on the McLerran-Venugopalan model [9], for the energy and momentum deposited in the central rapidity region by the colliding nuclei. Similar numerical results were obtained in Refs. [10][11][12]. This produced matter is oftentimes referred to as glasma; it is the precursor to quark-gluon plasma. We consider the backreaction of the glasma on the receding nuclei. Imposing energy and momentum conservation on the whole system, the increase of energy and momentum for the glasma is equal to the descrease of energy and momentum of the nuclei. Solving the energy momentum conservation equations gives us the anticipated rapidity loss and excitation energy for the receding nuclear fireballs. We neglect transverse coupling over the brief time duration of the collision, which is less than 1 fm/c. This means that the collisions are treated as a sum of independent tube-tube or streak-streak collisions, dependent on the transverse position r ⊥ . Once we have the baryon density and energy density, we can estimate the resulting temperature and baryon chemical potential of the fireballs as functions of r ⊥ . This is done with the help of a realistic crossover equation of state arXiv:1808.05751v1 [nucl-th] 17 Aug 2018 that smoothly connects the quark-gluon plasma phase at high temperature and high baryon chemical potential to the hadronic resonance gas phase at low temperature and small baryon chemical potential [13]. Apart from the prototype Au+Au collision at √ s N N = 200 GeV, we also explore other collision configurations, such as the asymmetric Cu+Au collision and collisions involving distorted nuclei like U+U collisions. Non-central collisions are also studied. In addition, we consider Pb+Pb collisions at even higher energy at the LHC. The outline of the paper is as follows. In Sec. II the formalism for implementing energymomentum conservation is presented. Properties of the equations of motion and determination of the input parameters are discussed. Section III is devoted to the prototype Au+Au central collision at √ s N N = 200 GeV. In Sec. IV, systematic studies of the high baryon densities in high energy heavy-ion collisions are given. We consider the nuclear size dependence, the collision energy dependence, and the impact parameter dependence. Conclusions are presented in Sec. V.

II. BARYON RAPIDITY LOSS AND NUCLEAR EXCITATION ENERGY
In this section we discuss the equations of motion and input parameters, followed by some numerical results.

A. Equations of Motion
For high energy heavy-ion collisions in the center-of-momentum frame, the nuclei are highly Lorentz contracted spheres (or ellipsoids if the ground state nucleus is deformed). We neglect transverse dynamics, which should not be important during the fraction of a fm/c time interval of relevance. Then the collision can be thought of as a sum of tube-tube or streak-streak collisions, each taking place at a particular value of the transverse coordinate r ⊥ with the beam along the z-axis. The projectile tube has a 4-momentum per unit area in the center-of-momentum frame denoted by P µ P = (E P , 0, 0, P P ). The tube loses energy and momentum to the classical color electric and magnetic fields produced in the region between the two receding tubes. This loss is quantified by dP µ P = −T µν glasma dΣ ν (1) where dΣ ν = (dz, 0, 0, −dt) is the infinitesimal four-vector perpendicular to the hypersurface spanned by dt, dz, and unit transverse area. The energy-momentum tensor of the glasma has been calculated in Refs. [7,8] as The A and B are functions of proper time τ = √ t 2 − z 2 (and other input parameters to be discussed in the following), while the dependence on space-time rapidity η = 1 2 ln[(t + z)/(t − z)] follows from the fact that T µν glasma is a second-rank tensor in a boost-invariant setting. The longitudinal position of the tube z P is a function of time, z P = z P (t). The z P is related to the time t via the velocity v P = dz P /dt = tanh y P , where y P is the momentum-space rapidity of the tube. Hence all the quantities solely depend on t. Note that T µν glasma must be evaluated on the trajectory of the tube. Explicitly dE P (t, z P ) = −T 00 glasma (t, z P )dz P + T 03 glasma (t, z P )dt , dP P (t, z P ) = −T 30 glasma (t, z P )dz P + T 33 glasma (t, z P )dt . (3) It is useful to define the Lorentz invariant effective mass per unit area M P via the relations E P = M P cosh y P and P P = M P sinh y P so that dE P = cosh y P dM P + M P sinh y P dy P , dP P = sinh y P dM P + M P cosh y P dy P .
We then express the differential form of the energy momentum conservation in Eqs. (3) in terms of the Milne coordinates (τ, x ⊥ , η). Using the transformations τ = t 2 − z 2 P and η P = 1 2 ln[(t + z P )/(t − z P )], the pseudorapidity of the projectile slab η P follows the equation where the auxilliary relations between (t, z P ) and (τ, η P ) have been used. Substituting Eqs. (2), (4) and (6) into Eq. (3), we obtain the equations of motion for the rapidity y P and the effective mass M P of the projectile slab in terms of Milne coordinates Note that Eq. (7) should be supplemented by the equation for pseudorapidity η P via Eq. (5).
It is important to point out that y P and M P are dynamical variables while η P is a geometric variable coming from the coordinate transformation. Equations (5) and (7) are coupled first order differential equations. Given the functions A(τ ) and B(τ ) from the glasma energy-momentum tensor, initial conditions are needed to solve for y P , M P and η P . But before discussing the initial conditions and numerical solutions to the equations, a few important properties of Eqs. (5) and (7) are worth noting.
• If B(τ ) = 0, the glasma energy-momentum tensor becomes diagonal, T µν glasma = diag(A, A, A, −A), which is like the string rope model [14,15]. In that case, dM P /dτ = 0 and the effective mass M P does not change with time. No energy is deposited in the nucleus. On the other hand, the rapidity y P decreases with time. As long as A(τ ) = 0, y P decreases until it becomes negative when the yo-yo type motion in string rope models appears. Therefore, a nonvanishing B(τ ) is necessary to incorporate nuclear excitation energy and prevent the appearance of yo-yo type motion. The glasma energy-momentum tensor predicts a nonvanishing B(τ ). Physically, the function A(τ ) represents contributions from the longitudinal chromo-electromagnetic fields E z and B z , while the function B(τ ) represents contributions from the transverse chromo-electromangetic fields E i and B i . Initially, at τ = 0, A = 0 and B = 0, so the transverse fields are zero while the longitudinal fields are nonzero. As the glasma evolves, B gradually increases while A gradually decreases until they become nearly equal; see Sect. IIc.
• If B(τ ) > 0, the effective mass M P (τ ) increases or decreases with time according to the sign of y P (τ ) − η P (τ ). For the excitation energy M P (τ ) to increase with time, y P (τ ) < η P (τ ) has to be maintained. From the definition tanh y P = dz P /dt, one obtains z P = t 0 tanh y P (t )dt > tanh y P (t) t because the velocity of the projectile tube v P = tanh y P decreases with time. With tanh η P = z P /t, the condition y P < η P is strictly guaranteed. Hence, positive values of B(τ ) predicts an increase of the effective mass and nuclear excitation energy. In addition, as long as y P = η P , the effective mass always increases with time, which brings up the question of a cutoff time dependence of the nuclear excitation energy. Physically, the increase of effective mass is governed by the strength of the transverse chromo-electromagnetic fields in the glasma which gradually builds up over time. It is worthwhile pointing out that this increase of transverse chromo-electromagnetic fields is not due to the decrease of the effective mass but comes from the decrease of the longitudinal chromo-electromagnetic fields.
• For y P < η P , Eq. (5) predicts that the pseudorapidity η P , like y P , also decreases with time.
This decrease of η P is purely kinematic in nature as a way to respond to any changes of y P . On the other hand, Eq. (7) predicts that the rate of momentum space rapidity loss dy P /dτ becomes smaller and smaller and finally saturates due to the time dependent properties of the functions A(τ ) and B(τ ).
To summarize, the nonvanishing, positive function B(τ ) in the glasma energy-momentum tensor, which is due to the transverse chromo-electromagnetic fields that gradually build up, is responsible for the nuclear excitation and prevents the momentum space rapidity from forever decreasing.

B. Determining the Input Parameters
Initial conditions must be specified. For nucleus-nucleus collisions at center-of-momentum collision energy per nucleon pair √ s N N , the initial rapidity y P (τ = 0) is computed by where m N = 0.937 (8) GeV is the nucleon mass. All tubes have the same initial momentum space rapidity. The initial value of pseudorapidity η P (τ = 0) is indeterminate, so we choose the initial value of η P to be equal to y 0 and start numerical calculations at an infinitesimal initial time τ = 0 + . Since we consider nucleus-nucleus collisions as a sum of tube-tube collisions at different r ⊥ , the initial effective mass, which is the mass per unit area, depends on the transverse position as Here T A (r ⊥ ) = dz ρ A (r ⊥ , z) is the nuclear thickness function for a nucleus with mass number A. The function ρ A (r ⊥ , z) is the Woods-Saxon distribution for a spherical nucleus.
Apart from the initial conditions, other input parameters come from the glasma energymomentum tensor. Specifically, the functions A(τ ) and B(τ ) in Eq. (2) factorize into two parts: the overall normalization ε 0 which is independent of time, and the time-dependent evolution functions. They can be expressed as The functions F A and F B depend on two free parameters, the ultraviolet cutoff scale Q and the infrared cutoff scale m on the transverse momentum. The ultraviolet cutoff scale characterizes the division between a description in terms of classical gluon fields and perturbative QCD . Larger values of Q attribute more energy and momentum to the classical fields whereas smaller values of Q attribute moere to the production of partons or minijets. The infrared cutoff is identified as the Λ QCD scale. The main complication comes from the initial energy density ε 0 . To facilitate further discussions, we quote its expression here from Ref. [8].
Clearly, for different slab-slab collisions characterized by different values of r ⊥ , the initial energy densities ε 0 (r ⊥ ) are different. This is due to the width µ(r ⊥ ) of color charge fluctuations per unit area. Here µ 1 (r ⊥ ) and µ 2 (r ⊥ ) are the color charge fluctuation widths for the two colliding slabs. Considering a nucleus, we assume that µ A (r ⊥ ) is a sum of contributions from all the nucleons at [16,17]. The µ N , which characterizes the gluon saturation for a nucleon, is related to the gluon saturation scale Q 2 s,N up to a logarithmic modification [18]. The proton saturation scale Q 2 s,N can be extracted from deep inelastic scattering experimental data [19]. With these considerations, we relate the color charge squared per unit area for a tube at transverse position r ⊥ to the tube at the central core region of the nucleus r ⊥ = 0 by As a consequence, the initial energy density for a tube-tube collision at transverse position r ⊥ scales as In arriving at Eq. (13), we assume the ultraviolet cutoff Q and the infrared cutoff m are the same for all the slabs given the same two colliding nuclei for a fixed collision energy. For asymmetric collisions, the respective nuclear thickness function has to be used. To determine the initial energy density at the central core ε 0 (r ⊥ = 0), initial conditions for hydrodynamic equations are invoked. Hydrodynamic equations require the initial value of the energy density ε hydro (τ 0 ) ≡ ε hydro (x = 0, y = 0; τ = τ 0 ) at spatial location x = 0, y = 0 with τ 0 the moment when hydrodynamics begins. Both ε hydro (τ 0 ) and τ 0 are free input parameters of hydrodynamic simulations and are tuned to reproduce the experimental data. We assume the initial classical gluon fields from the glasma state are valid until τ = τ 0 when it is switched to the hydrodynamic state. This is in the same spirit as the IP-Glasma model [16,17]. Therefore, the glasma energy density ε glasma = A(τ ) + B(τ ) (see Eq. (2)) at the central core equals ε hydro (τ 0 ) at τ = τ 0 .

C. Numerical Results
In this subsection we consider Au+Au collisions at zero impact parameter at the center-ofmomentum collision energy √ s N N = 200 GeV. We focus on the tube-tube collision that comes from the central core of the nucleus characterized by r ⊥ = 0 and compute the rapidity loss and nuclear excitation energy. The initial beam rapidity is y 0 = 5.36. The infrared cutoff is chosen to be m = Λ QCD = 0.2 GeV. The initial mass per unit area is M P (r ⊥ = 0) = 2.03 GeV/fm 2 . We use the hydrodynamic initial energy density ε hydro (τ 0 = 0.6 fm/c) = 30.0 GeV/fm 3 , which has been used in viscous hydrodynamic simulations in Ref. [20]. Depending on the ultraviolet cutoff chosen, ε 0 (r ⊥ ) has values 123. With these input parameters we solve the equations of motion Eqs. (7) and (5) up to proper time τ = 0.6 fm/c. Figure 3 shows the momentum-space rapidity y P of the central core of a gold nucleus as a function of proper time τ . The final rapidity in the center-of-momentum frame is 2.47. The central core loses about 2.9 units of rapidity within the first 0.1-0.2 fm/c; this is a robust result, insensitive to the value of Q. The rate of rapidity loss decreases with time at the beginning and finally approaches values close to zero, unlike the forever decreasing rapidity in string models. Figure 4 shows the excitation energy per baryon in units of the nucleon mass as a function of proper time. There is a slow but monotonic increase, unlike the rapidity loss whose asymptotic limit is reached within a few tenths of a fm/c. There is a weak dependence on Q. As can be seen from Eq. (7), the increasing rate of the nuclear excitation energy is determined by 2B(τ ) sinh (η P − y P ). At late time, B(τ ) maintains an almost constant postive value as shown in Fig. 2 while the difference η P − y P slowly diminishes but maintains finite positive value as shown in Fig. 5. As a consequence, the increasing rate of the excitation energy gradually decreases. At The ultraviolet cutoff is Q = 4.0 GeV. The condition y P < η P is maintained in the physically relevant proper time range.

III. HIGH BARYON AND ENERGY DENSITIES ACHIEVABLE AT TOP RHIC ENERGY
In the previous section we calculated the rapidity loss and excitation energy for one single tube-tube collision that comes from the central core of the nucleus. In this section, we repeat the calculations for all the tube-tube collisions at different r ⊥ with the help of Eqs. (9) and (13) for central Au+Au collisions at √ s N N = 200 GeV. Different tube-tube collisions characterized by r ⊥ will produce different rapidity losses y P (r ⊥ , τ ) and different excitation energies M P (r ⊥ , τ ). Once y P (r ⊥ , τ ) and M P (r ⊥ , τ ) have been obtained, the average rapidity loss is readily computed, which is constrained by the experimental data from the BRAHMS collaboration [21,22]. Furthermore, rapidity loss and excitation energy predict the baryon density and energy density that can be achieved. Additional thermodynamic properties of the high baryon density matter like temperature and baryon chemical potential can be calculated with the help of an appropriate equation of state.

A. Net-Baryon Rapidity Distribution
By collecting all the baryons at different r ⊥ after the collision at proper time τ 0 and record their final rapidity y P (r ⊥ , τ 0 ), one obtains the net-baryon rapidity distribution at that time. Let y P (r ⊥ ) denote the rapidity for the projectile tube at r ⊥ at proper time τ 0 = 0.6 fm/c. The netbaryon rapidity distribution dN B /dy is computed by summing up all the baryons at different r ⊥ that have rapidity y.
Here R A is the radius of the nucleus and and T A (r ⊥ ) the nuclear thickness function. We assume azimuthal symmetry in the transverse plane. The total number of baryons should be equal to the nuclear mass number A = y 0 0 dN B dy dy with y 0 the initial beam rapidity. The average rapidity loss follows from Figure 6 shows the net-baryon rapidity distribution at τ = 0.6 fm/c after the collision of gold nuclei at √ s N N = 200 GeV. The initial beam rapidities are y 0 = ±5.36 and the final rapidities for the central core of the Au nuclues are y P (r ⊥ = 0) = ±2.47. The central core of the gold nucleus experiences the largest rapidity loss while the peripheral part (r ⊥ ∼ R A ) experiences the smallest rapidity loss. For now we ignore possible thermal motion of baryons inside the nuclear tubes so that all the baryons at r ⊥ have the same rapidity y P (r ⊥ ). That is why there is a sharp vertical line at y = 2.47. The average rapidity loss is computed to be δy ≈ 2.4. The BRAHMS collaboration [21,22] was the only detector at the RHIC that could measure particle production anywhere near the fragmentaion regions. The coverage was limited to |y| ≤ 3.1, so the uncertainty in the rapidity loss estimate was large. For 0-5% centrality BRAHMS found an average rapidity loss of about 2.05 + 0.4/0.6. This is consistent with our result, especially since we focus on 0% centrality for illustration.

B. Large Baryon and Energy Densities
With the rapidity loss y P (r ⊥ , τ ) and the excitation energy M P (r ⊥ , τ ), we are able to calculate the baryon density and the energy density. The baryon density at proper time τ 0 is computed by [3][4][5][6] n B (r ⊥ , z , τ 0 ) = e ∆y(r ⊥ ,τ 0 ) ρ A (r ⊥ , z e ∆y(r ⊥ ,τ 0 ) ), where z = z − z P (r ⊥ , τ 0 ) is the longitudinal coordinate in the local rest frame of the tube characterized by r ⊥ . The change of rapidity ∆y(r ⊥ , τ 0 ) = y 0 − y P (r ⊥ , τ 0 ) depends on the transverse position r ⊥ . Different parts of the colliding nucleus characterized by r ⊥ have different time dependent rapidity y P (r ⊥ , τ ) and different longitudinal trajectory z P (r ⊥ , τ ) as viewed in the center-ofmomentum frame. Those tubes that are close to the central core of the nucleus experience larger rapidity loss and travel a smaller longitudinal distance z P (r ⊥ ∼ 0, τ 0 ) at time τ 0 , while those tubes close to the peripheral region of the nucleus experience small rapidity loss and travel to a larger longitudinal distance z P (r ⊥ ∼ R A , τ 0 ). Therefore, different tubes characterized by different r ⊥ are separated along the longitudinal direction due to different z P (r ⊥ , τ 0 ) at τ 0 , and the spherical shape of the nucleus before the collision will no longer be maintained after the collision. Since y P (r ⊥ , τ ) also depends on r ⊥ , there is no single reference frame that is the local rest frame for all the tubes comprising the fireball. Each tube has its own local rest frame by boosting the center-ofmomentum frame to the frame moving at rapidity y P (r ⊥ , τ 0 ). Multiplying the baryon density by the nuclear excitation energy, one obtains the energy density The energy density relies on the excitation energy which slowly increases with time, see Fig. 4. Hence, the energy density depends on the proper time chosen to evaluate its value. Figure 7 shows the proper energy and baryon densities as functions of the transverse coordinate at τ = 0.6 fm/c for z = 0. It should be noted that the maximum baryon density, about 3 baryons/fm 3 , is 20 times greater than the normal matter density of 0.155 nucleons/fm 3 . The maximum energy density is about 20 GeV/fm 3 , much larger than the critical energy density ∼ 1.0 GeV/fm 3 for the formation of quark-gluon plasma. Figure 8 is a contour plot of the proper baryon density. The contours are drawn at n B = 3, 2, 1, 0.5, and 0.15 baryons/fm 3 . The shapes of the contours arise for the following reasons. The diameter of a gold nucleus 2R A is about 14 fm. The core centered at r ⊥ = 0 along the z axis contains the most matter, suffers the greatest deceleration, and hence the greatest compression. Moving outward with increasing r ⊥ , the length of the tube is decreased to 2 R 2 A − r 2 ⊥ , and the deceleration, and hence compression, is reduced. These opposing effects approximately cancel each other, giving rise to roughly rectangular contours in the r ⊥ − z plane. Care must be taken when interpreting this figure. Since the rapidity loss depends on r ⊥ it means that there is a shear in the r ⊥ direction, and there is no single global frame of reference for all elements of the fireball. It should be emphasized that the baryon densities calculated here are more robust than the energy densities. The reason can be seen by comparing Eqs. 3 and 4. The rapidity loss, and therefore compression, is determined mostly within the first few tenths of a fm/c when the glasma dominates the dynamics. The excitation energy continues its slow growth as time goes on. If the transition from glasma to quark-gluon plasma happens earlier than 0.6 fm/c, it would reduce the excitation energy but hardly affect the compression. Exactly how the transition occurs is a topic of much current interest and activity. This should be kept in mind in the following discussions. The numbers are in units of baryons per fm 3 . The horizontal axis measures the distance along the beam direction in the local rest frame. Care must be taken when interpreting this plot since the rapidity of the matter and therefore the frame of reference depend on r ⊥ .
Both the baryon density Eq. (17) and the energy density Eq. (18) are evaluated in the local rest frame of the tube characterized by r ⊥ . In this representation, Fig. 8 displays the volume information of the fireball. To present the distributions in the center-of-momentum frame, where the spatial distribution of baryons is apparent, requires making Lorentz boosts from the individual local rest frame characterized by r ⊥ to the center-of-momentum frame. Note that the baryons shown in Fig. 8 share the same proper time τ instead of the same Cartesian coordinate time t. Figure 9 is a schematic picture of the central tube after collision in its own rest frame. Let (t * , z * ) denote Cartesian coordinates in the center-of-momentum frame and (t , z ) denote Cartesian coordinates in the local rest frame of the nuclear slab. Making the Lorentz transformation one gets where γ P = cosh y P and β P = tanh y P with y P the rapidity of the tube. The baryons all have the same proper time τ when viewed in the center-of-momentum frame so that z * L = τ sinh η L and t * L = τ cosh η L . Here η L is the pseudorapidity defined in the center-of-momentum frame. Then z L = τ sinh (η L − y P ) which, compared with z * = τ cosh η L , is just a shift of rapidity from η L to η L − y P . Therefore, the space-time pseudorapidity η L of a baryon labeled by L in the center-of- momentum frame is related the the corresponding coordinate z L in the local rest frame by A potential problem with Eq. (20) is that for z O = 0 it predicts η O = y P , which is not exactly true since η O slightly deviates from y P ; see Fig. 5. Here O indicates the center of the tube in Fig.  9. However, this only slightly influences the absolute position of the pseudorapidity for the center of the tube; the span of the pseudorapidity remains unchanged. This analysis improves upon and supercedes that reported in Ref. [4]. Figure 10 shows the baryon distribution in the r ⊥ − η plane using Eq. (20). The central tube spans about 1.5 units of rapidity. This distribution is useful as an initial condition for the subsequent hydrodynamic evolution in space and time, which is outside the scope of this paper.

C. Temperature and Baryon Chemical Potential
The previous results and discussions have not assumed the fireballs in the fragmentation regions to be thermalized. High energy densities are only necessary but not sufficient conditions for the formation of quark-gluon plasma. If and when thermalization occurs is hard to know. The proper time τ = 0.6 fm/c is the typical time when quark-gluon plasma is believed to be formed in the central region of Au+Au collisions at √ s N N = 200 GeV and is also the time when the matter in the central region is assumed to be thermalized. Exactly how the quark-gluon plasma in the central region equilibrates and thermalizes has not reached a definitive conclusion yet and is still under active investigation [23]. The practical approach is to assume the thermalization has been reached and tune the input parameters at τ = 0.6 fm/c (for example) so as to reproduce experimental data. The exact mechanism for the thermalization is actually not so relevant as long as it predicts the required initial conditions at τ = 0.6 fm/c constrained by experimental data. As for the fragmentation regions, Anishetty, Koehler and McLerran [3] argued that the matter in the fireballs in the fragmentation regions could thermalize due to interactions of secondary particles.
Since then, there has been very little research concerning the fragmentation regions, not to mention the thermalization problem. Thermalization in the fragmentation regions, if it is possible, could be as challenging as the thermalization problem in the central region and beyond the scope of this paper. Just like the current practice in the central region, we assume the fireballs have reached thermalization and explore the consequences. Assuming local thermal equilibrium, what are the temperature and baryon chemical potential? To answer that question we need an equation of state. In the limit of very high temperature and baryon density, the equation of state can be calculated using perturbative QCD in thermal field theory [24,25]. On the other hand, in the low temperature and small baryon density regime, the relevant degrees of freedom are hadrons. For a thermalized system of hadrons, the hadron resonance gas model gives a very good description of the thermodynamic properties of the system. Furthermore, first principle calculations based on lattice QCD provides robust results for the equation of state of a system of quarks and gluons in a very wide range of temperature for zero baryon chemical potential [26,27]. Extending to finite baryon chemical potential has the notorious sign problem. Currently, much effort has been devoted to extending the lattice calculations to finite baryon chemical potential [28]. In the following, we will use a crossover equation of state [13] that smoothly connects the quark-gluon plasma phase and the hadronic resonance gas phase consistent with lattice data [27,29]. This equation of state does not contain a first order phase transition line or critical point. Instead, the transition from the quark-gluon plasma phase to the hadronic resonance gas phase is a rapid, smooth crossover both for zero baryon chemical potential and for nonzero baryon chemical potentials. The crossover equation of state [13] has the form with the switching function Here P qg represents the perturbative QCD results for the pressure of the quark-gluon plasma phase while P h represents pressure from the excluded volume model of the hadronic resonance gas phase. The switching function S(T, µ) asymptotically approaches 1 for very large T and µ and asymptotically approaches 0 for very small T and µ. The free parameters T 0 , µ 0 and r are optimized to be consistent with lattice data. Given the energy and baryon density distributions for z = 0 shown in Fig. 7, we compute the corresponding temperature and baryon chemical potential distributions using the crossover equation of state from Eqs. (21) and (22). Figure 11 shows the temperatures and baryon chemical potentials for different rapidities y P instead of r ⊥ . Note that there is a one-to-one correspondence between r ⊥ and y P : y P = y P (r ⊥ ). Only those areas with energy densities larger than 1.0 GeV/fm 3 have been displayed. In the range of rapidity from 2.47 to 3.45, the value of the baryon chemical potential increases from 655 MeV to 1020 MeV while the value of the temperature correspondingly decreases from 328 MeV to 155 MeV. The baryon chemical potential to temperature ratio ranges from 2.0 to 6.5 as shown in Fig. 12. The same data is shown in the µ B − T plane in Fig. 13. Assuming the subsequent expansion of the thermalized high baryon density matter after τ = 0.6 fm/c is approximately adiabatic, just like in the central region, then typical phase trajectories in the µ B − T plane may be computed. See Fig. 14  ratios. The corresponding momentum space rapidities are correlated with the entropy per baryon. These phase trajectories tilt to the right after changing from the quark-gluon plasma phase to the hadronic gas phase. In the transition region, the temperature decreases very quickly while the baryon chemical potential decreases very slowly. From these phase trajectories, it is possible that the expansion of the high baryon density matter in the fragmentation regions might go through or near the region in the µ B − T plane where a first order phase transition line or a critical point are conjectured to occur [30,31]. As shown in Fig. 15, the entropy per baryon ratio might be in the right range so that a scan through the momentum space rapidity may locate the critical point.

IV. HIGH BARYON DENSITY ACHIEVABLE IN OTHER COLLISION CONFIGURATIONS
To extend the previous analysis of the high baryon densities achieved outside the central rapidity region of Au+Au collisions at √ s N N = 200 GeV to other heavy ion collision configurations involving different nuclear sizes, different collision energies or different impact parameters, the crucial step is to determine the initial energy density of the glasma ε 0 (r ⊥ ) for these different collision configurations. In the case of central Au+Au collisions at √ s N N = 200 GeV, the ε 0 (r ⊥ ) is determined by ε 0 (r ⊥ = 0) through the nuclear thickness function T A (r ⊥ ), while ε 0 (r ⊥ = 0) is fixed by the initial energy density of hydrodynamics at τ 0 = 0.6 fm/c, assuming that the classical gluon fields are valid up to the time when hydrodynamic evolution of the quark-gluon plasma begins. In principle, we can repeat this procedure to determine ε 0 (r ⊥ ) for other collision configurations as long as we know the starting time of hydrodynamics τ 0 and the energy density at that time ε(r ⊥ = 0, τ 0 ). These two quantities are usually optimized to reproduce bulk experimental data after running hydrodynamic simulations. In the following analysis, however, we take a different approach in determining ε 0 (r ⊥ ) for other collision configurations. We assume that Eq. (11) is the formula which determines the initial energy density ε 0 (r ⊥ ). Rewritten more explicitly, the initial energy density is a functional of µ i (r ⊥ , √ s) and Q i (r ⊥ , √ s) with i = 1, 2 indicating the two colliding nuclei. We assume the infrared cutoffs to be the same m 1 = m 2 = Λ QCD for the two colliding nuclei. Both the color charge squared per unit area µ i and the ultraviolet cutoff Q i depend, in principle, on the transverse position r ⊥ and the center-of-momentum collision energy √ s. The color charge squared per unit area for a nucleus is related to that of a nucleon by We assume, for a given nucleus-nucleus collision, that the ultraviolet cutoff is independent of the transverse position r ⊥ so that The central Au+Au collision at √ s N N = 200 GeV will serve as a reference for other collisions involving different nuclear sizes, different collision energies and different impact parameters.

A. Nuclear Size Dependence
In this subsection, we consider heavy-ion collisions at fixed center-of-momentum collision energy √ s N N = 200 GeV. We will study Cu+Cu, Cu+Au, and U+U collisions in addition to Au+Au collisions. RHIC has already run Cu+Cu, Au+Au, and Cu+Au collisions at √ s N N = 200 GeV and run U+U collisions at √ s N N = 193 GeV. For a fixed center-of-momentum collision energy, the initial energy density for a general nucleus-nucleus central collision can be obtained with reference to the Au+Au collision by The initial energy density for Au+Au collisions is Computation of nuclear thickness functions involves the Woods-Saxon nuclear distribution function. For spherical nuclei, such as Au and Cu, the Woods-Saxon distributions ρ A (r) are spherically symmetric. For deformed nuclei like U, we use the following parameterization [32,33] Here ρ 0 = 0.166 fm −3 , ξ = 0.44 fm, R A = 6.86 fm and Ω(θ) = 1 + β 2 Y 2 0 (θ) + β 4 Y 4 0 (θ) with β 2 = 0.280 and β 4 = 0.093. The Y l m (θ) are spherical harmonic functions. The angle θ is related to the Cartesian coordinates by sin θ = r ⊥ /r and cos θ = z/r. The θ = 0 corresponds to the direction of the longest axis while θ = π/2 corresponds to the direction of the shortest axis. Equation (27) describes an ellipsoid-like shape. In central U+U collisions, depending on the orientations of the uranium nuclei, there could be many different collision configurations. In the following discussion, we only consider the tip-tip collision configuration where the longest axes of the uranium nuclei align with the beam directions.  Here ε hydro is the initial energy density in the central rapidity region with x = y = 0 at τ = 0.6 fm/c when hydrodynamics starts. The average rapidity loss is δy . The y P denotes the final rapidity of the central core of the nucleus, which experiences the largest rapidity loss. Other thermodynamic quantities n B , ε P , T , µ B and s are also evaluated in the central core of the fireball.
all have the same initial beam rapidity y 0 = 5.36. The average rapidity loss δy increases with the increase of nuclear mass. This can be understood from the tube-tube collision at the central cores of the colliding nuclei where the remaining quantities y P , n B , ε P , T and µ B are evaluated.
In the three collisions, the central core of a projectile starts with initial rapidity y 0 and ends with rapidity y P . For collisions involving nuclei of larger atomic mass, like U+U, the final rapidity y P is smaller and the rapidity loss experienced is therefore larger. This point is encoded in the equations of motion (7) where M P ∼ T A (r ⊥ ) and A, B ∼ [T A (r ⊥ )] 2 so that the rate of rapidity change dy/dτ ∼ T A (r ⊥ ). When averaging over r ⊥ , collisions of nuclei with bigger nuclear size have larger average rapidity loss. The average rapidity loss in Au+Au collisions at √ s N N = 200 GeV for the centrality class 0 − 5% has been measured and estimated to be in the range from 1.45 to 2.45 by the BRAHMS collaboration [21]. However, the average rapidity losses for the Cu+Cu collision and the U+U (tip-tip) collision have not been measured experimentally. Baryon densities increase with increasing nuclear mass in accordance with the rapidity losses because the baryon density is proprotional to the exponential of the rapidity loss. These maximal baryon densities are all more than ten times larger than the normal nuclear density. The maximal energy densities ε P obtained in the fragmentation regions of the three collisions are smaller than the respective energy densities ε hydro in the central rapidity region at τ = 0.6 fm/c when hydrodynamics begins. This is consistent with the expectation that energy density achieved in the central region of high energy heavy-ion collisions is larger than the energy density achieved in the receding nuclear fireballs. The temperature T and baryon chemical potential µ B corresponding the largest baryon density n B and energy density ε P for the three collisions are also given in Table I. Finally, the largest entropy per baryon s/n B that can be achieved in the fragmentation regions of the three collisions are 18.33, 23.53 and 27.82 which increases with the nuclear atomic mass, reflecting the fact that the energy density increases faster than the baryon density.  Figures 16, 17, 18 and 19 show the baryon density, energy density, baryon chemical potential and temperature achievable in the central core of the three collisions as functions of r ⊥ . The temperature and baryon chemical potential are displayed only in the regions with energy density larger than 1.0 GeV/fm 3 . In Fig. 18, the baryon chemical potential increases as r ⊥ increases, reaches a maximum, and then decreaes. This feature is due to the equation of state because the corresponding temperatures in these regions are in the transition regions from the quark-gluon plasma to the hadronic resonance gas. Figures 20 and 21 show the baryon density distributions in the center-of-momentum frame. Here η is the space-time pseudorapidity. For Cu+Cu collisions, the pseudorapidity spans about 1.5 unit around the central core regions while for U+U collision, the pseudorapidity spans about 1.0 unit. For asymmetric central collisions of Cu+Au at √ s N N = 200 GeV, the gold nucleus has a larger diameter than the copper nucleus so that the transverse overlap size is only the cross sectional area of the copper nucleus. Peripheral regions of the gold nucleus play the role of spectators and do not participate in the collision. To solve the equations of motion (7) for the copper and gold fireballs after the collision, one has to generalize the initial energy density in Eq. (26) to incorporate contributions from two different nuclei as Here A and A indicate the two colliding nuclei. As a consequence, the functions A, B in the equations of motion (7) scale as T A T A . For the copper nucleus, after canceling T Cu from the two sides of equations in (7), the equations of motion only depend on the thickness function T Au . These equations of motion are exactly the same as the equations of motion for the gold nucleus in symmetric Au+Au collisions at the same energy. Likewise, the equations of motion for a gold nucleus in Cu+Au collisions are exactly the same as the equations of motion governing the copper nucleus in Cu+Cu collision at the same energy. In other words, concerning the rapidity loss and the excitation energy, the following equivalences are valid.  GeV; this is the initial distribution to be used in hydrodyanic modeling. Peripheral regions of the gold nucleus do not participate in the collision so that their rapidities are the same as the beam rapidity y 0 = −5.36. The net baryon distribution on the copper side is smaller magnitude and narrower when compared to the net baryon distribution in Au+Au collisions as shown in Fig.  6, but they have the same smallest rapidity cutoff at y = 2.47.

B. Collision Energy Dependence
In this section we study the high baryon densities achievable with different collision energies. Specifically, we focus on Au+Au central collisions at fields becomes more applicable with increasing center-of-momentum collision energy. For the Au+Au collision at √ s N N = 62.4 GeV, whether the McLerran-Venugopalan model is still applicable needs more detailed study, which is beyond the scope of this paper. Here we explore the Au+Au collision at √ s N N = 62.4 GeV as an extrapolation of our method. Importantly, the average baryon rapidity loss and net proton rapidity distribution have been measured for Au+Au collision at √ s N N = 62.4 GeV [22].
To solve for the rapidity loss and excitation energy from Eq. (7) for collisions at different energies, one needs to generalize the calculation of initial energy density from Eq. (26) to include the collision energy dependence. With the help of Eq. (24) we have In obtaining this expression, we ignore the expected collision energy dependence of the ultraviolet cutoff illustrated in Eq. (25) which appears as ln Q 2 /m 2 in the initial energy density expression of Eq. (11). In principle the ultraviolet cutoff Q has to be adjusted when the saturation scale Q s changes with collision energy. This is to make sure that the scale separation Q 2 Q 2 s is satisfied so that the semi-analytic expression of the glasma energy-momentum tensor we obtained in the leading Q 2 approximation can be used. In Au+Au collisions at √ s N N = 200 GeV, we used the value Q s = 1.2 GeV and Q = 4.0 GeV. As it will become clear in the following, the saturation scale at the LHC energy 5.02 TeV is approximately Q s ∼ 1.8 GeV and an ultraviolet cutoff of Q ∼ 6.0 GeV only contributes approximately 10% change in the energy density after taking the logarithm. This logarithmic change of the ultraviolet cutoff in the expression for the energy density should be minor compared to the power law changes in the µ N ( √ s N N ). Note that µ N is related to the saturation scale of the nucleon Q 2 sN up to a logarithmic correction. The collision energy dependence of the saturation scale Q 2 sN (x) can be parameterized as [34][35][36] with Q 2 0 = 1.0 GeV 2 , x 0 = 3 × 10 −4 and λ = 0.288. This parametrization provides a successful description of the HERA data on deep inelastic scattering. At central rapidity y = 0, the small x parameter is related to the center-of-momentum collision energy by x = Q sN / √ s N N , and one obtains [37] Q 2 sN ( √ s N N ) whereλ = λ 1+λ/2 = 0.252 with √ s 0 some reference collision energy. Therefore, the final expression for the collision energy dependent initial energy density is Using this simple parametrization of the collision energy dependence, we compute the initial energy density of the glasma and the hydrodynamic initial energy density for the different collision energies; the results are given in Table III.

C. Non-Central Collisions
Apart from colliding heavy ions of different sizes and collision energies, realistic heavy-ion collisions also measure observables at different centralities. The centrality characterizes the impact paramter of the two colliding nuclei. From the most central collisions with impact parameter b ∼ 0 to the most peripheral collisions with impact parameter b ∼ 2R A , the centrality increases from 0% to 100%. The impact parameter, however, cannot be directly measured in heavy-ion experiments. It is to be noted that the observed particle multiplicity decreases from central collisions to peripheral collisions. Therefore, experiments measure the particle multiplicity and categorize these events using centrality values in the range 0% − 100%. With the help of a geometric model, like the Glauber model, one can then relate the centrality to the impact paramter. A few examples are given in Table V as calculated by the PHENIX collaboration [39]. In this subsection, we will explore the impact parameter dependence of the baryon densities in the fireballs in high energy heavy-ion collisions.  To characterize non-central collisions in the transverse plane we use cylindrical coordinates. Figure 26 is a schematic illustration of the transverse overlap region of two equal size nuclei colliding at non-zero impact paramter b. All the vectors in this figure are 2-dimensional in the x-y plane. Let the projectile P be located at b/2 and the target T at −b/2 with b = (0, b). An arbitrary point is labeled by r ⊥ = r ⊥ (cos φ, sin φ) = (x, y). The distance from the centers of the nuclei to that point are The thickness functions T P and T T depend only on r P and r T , respectively. Consider the baryon distribution arising from the projectile; a similar formula applies to the target. Let y P (r ⊥ , φ, b) denote the final rapidity of a piece of projectile located at the position r ⊥ . Then the baryon rapidity distribution is a generalization of Eq. (15) taking into account the non-zero impact parameter, For central collisions of spherical nuclei, which need not be identical, there is no dependence on φ, and there is a one-to-one correspondence between y P and r. Then we can replace δ(y − y P (r ⊥ )) with δ(r ⊥ − r ⊥P (y)) along with the relevant Jacobian to reproduce the expression in Eq. (15). For b = 0 it is may be better to discretize both r and φ. Assume uniform spacings ∆r and ∆φ which satisfy N r × r =R and (N φ + 1) × ∆φ = 2π, respectively. HereR characterizes the radius of a circular area chosen for the discretization. The discretized points are labeled by (r i , φ j ) with r i = i∆r, i = 1, 2, . . . , N r φ j = j∆φ, j = 0, 1, 2, . . . , N φ . The point (0, 0) is treated separately. For points in the overlap region, their coordinates (r i , φ j ) satisfy r 2 P(i,j) ≤ R 2 A and r 2 T(i,j) ≤ R 2 A . For the tube-tube collisions at the position (r i , φ j ), the mass per unit area is obtained by the thickness functions M P(i,j) = m N T A (r P(i,j) ) and M T(i,j) = m N T A (r T(i,j) ), while the initial energy density is proportional to the product of the two nuclear thickness functions ε 0(i,j) ∼ T A (r P(i,j) )T A (r T(i,j) ). We can then solve for the final rapidity y P(i,j) and the nuclear excitation energy M P(i,j) . The number of baryons within the projectile tube at position (r i , φ j ) is ∆N P(i,j) = r i T P(i,j) ∆r∆φ. To calculate the baryon rapidity distribution for b = 0, we perform an intergration over the Dirac delta function while fixing the value of φ.
In Figs. 27 and 28 the proper baryon density distribution of the fireball in the transverse plane for the longitudinal slice z = 0 is given for b = 0 and b = 2.3 fm, respectively. For zero impact parameter, the baryon distribution is azimuthally symmetric in the transverse plane and the center of the fireball achieves the largest baryon density of about 3 baryons/fm 3 . In contrast, for collisions with non-zero impact paramter b = 2.3 fm, the baryon distribution in the transverse plane is no longer azimuthally symmetric. Even the region where the largest baryon density is achieved has been shifted away from the center of the projectile fireball (the center of the fireball corresponds to (x = 0, y = 0)). The region with the largest baryon density surrounds the point (x = 0, y = 2.3 fm), which is the transverse location of the center of the target fireball when the projectile nucleus and the target nucleus overlap. This is easy to understand from our previous discussion on asymmetric Cu+Au collisions. When two nuclear tubes collide, the tube with less baryon charge experiences larger nuclear compression, while the tube with more baryon charge experiences less nuclear compression.  Figure 29 shows the baryon density distribution in the x − η plane, where the impact parameter is in the y direction. The region of the projectile nucleus that overlaps with the central region of the target nucleus experiences the largest nuclear compression while the regions of the projectile nucleus that overlap with the peripheral region of the target nucleus experiences less nuclear compression.

V. CONCLUSION
In this paper we have systematically studied the high baryon densities outside the central rapidity region of high energy heavy-ion collisions within the McLerran-Venugopalan model. The off-diagonal term in the energy-momentum tensor of the glasma, which comes from the transverse chromo-electromagnetic fields, is reponsible for the nuclear excitation energy in this model. In contrast, typical string models do not deposit energy in the receding nuclei. For central Au+Au collisions at √ s N N = 200 GeV, the highest baryon density 3.0 baryons/fm 3 is about twenty times larger than normal nuclear density, and the largest energy density 20 GeV/fm 3 is more than one hundred times larger than the energy density of nuclear matter at its saturation density. Using a crossover equation of state, the temperature ranges from about 155 to 330 MeV, while and the baryon chemical potential ranges from about 650 to 1020 MeV. The entropy per baryon ratios are found to be in the range of 10 to 23.5, corresponding to the momentum space rapidity range 2.5 to 3.7. The entropy per baryon ratios might be in the right range so that a scan through rapidity may locate the critical point of QCD phase diagram for central collisions. For central collision at fixed beam energy, the highest baryon density achievable in the fragmentation regions increases with nuclear size. In Cu+Cu, Au+Au and U+U (tip-tip) collisions at √ s N N = 200 GeV, the highest baryon densities are about 2, 3 and 4 baryons/fm 3 , respectively. In asymmetric Cu+Au collisions, the Cu nucleus is compressed more than the Au nucleus so that the Cu fireball achieves the higher baryon densities. For central collisions at fixed nuclear size, the highest baryon density achievable increases with collision energy. We numerically studied Au+Au collisions at √ s N N =62. 4

and 200
GeV and Pb+Pb collisions at √ s N N = 2.76 and 5.02 TeV. For non-central collisions, the average rapidity loss is reduced due to the spectator baryons in the peripheral region of the colliding nculei. For Au+Au collisions at √ s N N = 200 GeV with 0 − 5% centrality, our calculations predict the average rapidity loss to be about 2.14. This is in agreement with measurement by the BRAHMS collaboration. Rapidity loss at LHC is not known because of the challenge of measuring and identifying particles with large rapidity in the detector's frame of rest. We must emphasize that the results obtained here provide the initial conditions for relativistic fluid dynamic descriptions of high energy heavy ion collisions and cannot be compared directly to experimental data. Improvements can be made within the theoretical framework described here. For example, the McLerran-Venugopalan model for the glasma, as implemented by us, could be complimented by the inclusion of the production of minijets. An equation of state could be used that incoporates a critical point. But perhaps the biggest challenge is how to identify and measure baryons (and mesons) at the higher rapidities of relevance to the high baryon density matter.