3. Before hydrodynamics: initial conditions, energy deposition & pre-equilibrium dynamics¶
In these notes we define the early stage of heavy ion collisions as everything that happens before the plasma can be described by hydrodynamics.
Hydrodynamics describes the evolution of the strongly-coupled quark-gluon plasma in space and time. It can describe this evolution if it is provided the state of the fluid on a given spacetime surface [*].
| [*] | That is, if hydrodynamics is taken as nothing more than a system of partial differential equations, boundary conditions or initial conditions must be provided before a numerical solution of the equations can be obtained. |
Recall that, with second-order viscous hydrodynamics, the energy-momentum tensor \(T^{\mu\nu}(X)\) is decomposed into four independent hydrodynamic fields:
- The energy density \(\epsilon(X)\)
- The flow velocity \(u^\mu(X)\)
- The shear stress tensor \(\pi^{\mu\nu}(X)\)
- The bulk pressure \(\Pi(X)\)
In almost every simulation of heavy ion collisions, the equations of hydrodynamics are initialized at a fixed value of \(\tau\): \(\tau_0\). That is, hydrodynamic initial conditions are provided as \(T^{\mu\nu}(\tau=\tau_0,x,y,\eta_s)\). More typically the different fields are initialized separately: \(\epsilon(\tau=\tau_0,x,y,\eta_s)\), \(u^\mu(\tau=\tau_0,x,y,\eta_s)\), \(\pi^{\mu\nu}(\tau=\tau_0,x,y,\eta_s)\) and \(\Pi(\tau=\tau_0,x,y,\eta_s)\)separately.
Note that these initial conditions represent a considerable amount of information. To understand how much of a challenge it is to determine \(T^{\mu\nu}(\tau=\tau_0,x,y,\eta_s)\), it is worth taking a step back and look at the initial moments of a heavy ion collisions.
3.1. The earliest stages of heavy ion collisions¶
No two collisions of nuclei are the same, and one of the reasons can be understood from a simple geometric picture, illustrated in the following figure:
The two nuclei (ions) can overlap completely or can have a small overlap, and the strongly-coupled quark-gluon plasma that is formed in these two cases will be very different. The geometric degree of overlap can be quantified with the impact parameter “b”, shown on the above figure.
In general, the more the overlap (the more “central” the collision is, “b” small), the larger the plasma. Smaller overlap between the nuclei results in smaller and shorter-lived plasma, and while very peripheral heavy ion collisions that barely graze each other might not produce any quark-gluon plasma at all.
Note
The geometric picture discussed above is a very intuitive way of picturing heavy ion collisions. However, it is important to remember that the impact parameter “b” is not a quantity that can be extracted from heavy ion collision data. In this sense, it can be preferable to think of it as a model parameter, not a measurable quantity.
There are additional fluctuations aside from fluctuations in the degree of overlap of the nuclei: the plasma produced by two collisions with the same overlap will differ, since collisions of nuclei are quantum processes with a probabilistic rather than deterministic outcome.
The other major factor that affects the size and the lifetime of the plasma is the center-of-mass energy, \(\sqrt{s_{NN}}\). The higher the collision energy, the larger and longer-lived the produced quark-gluon plasma.
Note
MUSIC and hydrodynamic models in general have been found to describe a large number of experimental measurements in collisions of gold nuclei (Au) at the top RHIC energy (\(\sqrt{s_{NN}}=200~\textrm{GeV}\)) and collision of lead ions (Pb) at the LHC (\(\sqrt{s_{NN}}=2760~\textrm{, } 5020~\textrm{GeV}\), …). It is likely that hydrodynamic models can still be used in collisions of smaller nuclei and/or lower \(\sqrt{s_{NN}}\) collisions, but one should proceed with additional caution as these systems can have additional complications.
3.2. Obtaining initial conditions for hydrodynamics¶
From a purely practical point of view, we can ask the question
“Given the species of the colliding nuclei (e.g. Au or Pb) and a collision energy (\(\sqrt{s_{NN}}\)), what are the initial conditions \(T^{\mu\nu}(\tau=\tau_0,x,y,\eta_s)\) of hydrodynamics?”
This is a pragmatic question that hides a lot of the complexity of the early stage of heavy ion collisions. Between the moment that the nuclei collide and the time at which hydrodynamics become applicable, there is a complex dynamical evolution of the deconfined matter.
Moreover the fact that there are fluctuations in heavy ion collisions implies that there is no single initial condition corresponding to two nuclei colliding. There is rather a family, or a distribution, of initial conditions.
Describing this pre-hydrodynamic stage of stage of heavy ion is out of the scope of MUSIC. Broadly speaking, there are two options to initialise the hydrodynamics in MUSIC at \(\tau=\tau_0\):
- An external model can be used to describe the early stage of to collision. The results of this external model can be read in and used as initial conditions for MUSIC. This is the procedure adopted when IP-Glasma initial conditions are used with MUSIC.
- The initial conditions of hydrodynamics can be parametrized with the Glauber ansatz
The initialization of the hydrodynamics fields is performed in function “Init::InitTJb()” in file “init.cpp”. This function first allocates memory for the discretized hydrodynamic fields and then initialize the fields. The input parameter Initial_profile is used to choose which type of initial conditions are used to the subsequent hydrodynamic evolution. The most important options are
###################################
# parameters for initial conditions
###################################
Initial_profile 3 # type of initial condition
# 1: Optical Glauber model
# 3: MC-Glauber model
# 8: Read in initial profile from a file
# (e.g., IP-Glasma)
Reading-in an external initial condition is a more flexible option, but in certain cases Glauber initial conditions are sufficient.
3.2.1. Reading-in external hydrodynamic initial conditions¶
Since there is no accepted universal format of hydrodynamic initial conditions, reading-in external initial conditions will almost certainly require modifications to the code. This should be straightforward, however, given examples already in MUSIC that show how to initialize hydrodynamics with external files.
When the input parameter Initial_profile 8 is used, MUSIC will read an external file
specified by another input parameter, Initial_Distribution_Filename.
The format expected when using Initial_profile 8 is the following. First, there is a header which
indicates the number of grid points in each direction and the size of the grid cells, in fermi for x/y (unitless in \(\eta\)):
string dummy;
double dummy2;
profile >> dummy >> dummy >> dummy2
>> dummy >> neta >> dummy >> nx >> dummy >> ny
>> dummy >> deta >> dummy >> dx >> dummy >> dy ;
where “neta”, “nx” and “ny” are the number of grid points in each directions, and “deta”, “dx” and “dy” are the step size.
An example of header line would be:
# dummy 1 etamax= 1 xmax= 200 ymax= 200 deta= 0 dx= 0.17 dy= 0.17
The rest of the file is read as follows:
for (ix = 0; ix <= DATA->nx; ix++) {
for (iy = 0; iy <= DATA->ny; iy++) {
profile >> dummy1 >> dummy2 >> dummy3
>> density >> utau >> ux >> uy
>> dummy >> dummy >> dummy >> dummy;
where “density” is the energy density in GeV/fm^3, utau is u^tau, ux is u^x and uy is u^y. The “dummy” variable contain either redundant or unused information.
The information read from the external file is then used to initialize the discretized hydrodynamic fields:
int Init::InitTJb([...]) {
[...]
if (DATA->Initial_profile == 8) {
// read in the profile from file
// - IPGlasma initial conditions with initial flow
[...]
music_message.info(" ----- information on initial distribution -----");
music_message << "file name used: " << DATA->initName;
music_message.flush("info");
ifstream profile(DATA->initName.c_str());
string dummy;
int nx, ny, neta;
double dx, dy, deta;
// read the information line
profile >> dummy >> dummy >> dummy >> dummy >> neta
>> dummy >> nx >> dummy >> ny
>> dummy >> deta >> dummy >> dx >> dummy >> dy;
music_message << "neta=" << DATA->neta << ", nx=" << nx
<< ", ny=" << ny << ", deta=" << DATA->delta_eta
<< ", dx=" << dx << ", dy=" << dy;
music_message.flush("info");
[Allocates memory]
//read the one slice
for (ix = 0; ix <= DATA->nx; ix++) {
for (iy = 0; iy <= DATA->ny; iy++) {
profile >> dummy1 >> dummy2 >> dummy3
>> density >> utau >> ux >> uy
>> dummy >> dummy >> dummy >> dummy;
temp_profile_ed[ix][iy] = density;
temp_profile_utau[ix][iy] = utau;
temp_profile_ux[ix][iy] = ux;
temp_profile_uy[ix][iy] = uy;
if (ix == 0 && iy == 0) {
DATA->x_size = -dummy2*2;
DATA->y_size = -dummy3*2;
music_message << "eta_size=" << DATA->eta_size
<< ", x_size=" << DATA->x_size
<< ", y_size=" << DATA->y_size;
music_message.flush("info");
}
}
}
profile.close();
[Allocates more memory]
[...]
for(int ieta = 0; ieta < DATA->neta; ieta++)
{
double eta = (DATA->delta_eta*(ieta + DATA->neta*rank)
- (DATA->eta_size)/2.0);
double eta_envelop_ed = eta_profile_normalisation(DATA, eta);
for(ix = 0; ix <= DATA->nx; ix++)
{
for(iy = 0; iy<= DATA->ny; iy++)
{
rhob = 0.0;
[...]
epsilon = (temp_profile_ed[ix][iy]*eta_envelop_ed
*DATA->sFactor/hbarc); // 1/fm^4
[...]
if (epsilon<0.00000000001)
epsilon = 0.00000000001;
// initial pressure distribution
p = eos->get_pressure(epsilon, rhob);
// set all values in the grid element:
(*arena)[ix][iy][ieta].epsilon = epsilon;
(*arena)[ix][iy][ieta].epsilon_t = epsilon;
(*arena)[ix][iy][ieta].prev_epsilon = epsilon;
(*arena)[ix][iy][ieta].rhob = rhob;
(*arena)[ix][iy][ieta].rhob_t = rhob;
(*arena)[ix][iy][ieta].prev_rhob = rhob;
(*arena)[ix][iy][ieta].p = p;
(*arena)[ix][iy][ieta].p_t = p;
(*arena)[ix][iy][ieta].T = eos->get_temperature(epsilon, rhob);
(*arena)[ix][iy][ieta].mu = eos->get_mu(epsilon, rhob);
[Allocate even more memory]
/* for HIC */
u[0] = (*arena)[ix][iy][ieta].u[0][0] = temp_profile_utau[ix][iy];
u[1] = (*arena)[ix][iy][ieta].u[0][1] = temp_profile_ux[ix][iy];
u[2] = (*arena)[ix][iy][ieta].u[0][2] = temp_profile_uy[ix][iy];
u[3] = (*arena)[ix][iy][ieta].u[0][3] = 0.0;
(*arena)[ix][iy][ieta].prev_u[0][0] = temp_profile_utau[ix][iy];
(*arena)[ix][iy][ieta].prev_u[0][1] = temp_profile_ux[ix][iy];
(*arena)[ix][iy][ieta].prev_u[0][2] = temp_profile_uy[ix][iy];
(*arena)[ix][iy][ieta].prev_u[0][3] = 0.0;
(*arena)[ix][iy][ieta].pi_b[0] = 0.0;
for(int mu=0; mu<4; mu++)
{
/* baryon density */
(*arena)[ix][iy][ieta].TJb[0][4][mu] = rhob*u[mu];
for(nu=0; nu<4; nu++)
{
(*arena)[ix][iy][ieta].TJb[0][nu][mu] = (
(epsilon + p)*u[mu]*u[nu]
+ p*(DATA->gmunu)[mu][nu]);
(*arena)[ix][iy][ieta].Wmunu[0][nu][mu] = 0.0;
(*arena)[ix][iy][ieta].prevWmunu[0][nu][mu] = 0.0;
(*arena)[ix][iy][ieta].Pimunu[0][nu][mu] = 0.0;
(*arena)[ix][iy][ieta].prevPimunu[0][nu][mu] = 0.0;
}/* nu */
}/* mu */
}
}
}/* ix, iy, ieta */
// clean up
[Free up memory]
}
music_message.info("initial distribution done.");
return 1;
}
To summarize, there are two main input parameters used in this case:
Initial_profile 8 # 8: Read in initial profile from a file
# (e.g., IP-Glasma)
# external file from which the initial conditions are read
Initial_Distribution_input_filename example_inputs/filename
Extra parameters
For the example above, there are two extra parameters:
initialize_with_entropy 0 # 0: with energy density
# 1: with entropy density
s_factor 1.0 # normalization factor read in
# initial data file
The code in which parameter initialize_with_entropy appears was hidden from the previous example (but can be seen in the original source code), as I think this parameter should not need to be used when reading external files.
Parameter s_factor is designed to be used to normalize the energy density read-in from the external file Initial_Distribution_input_filename. Ideally it is not a parameter that should need to be used and it should be set to “1.0”.
3.2.2. Glauber initial conditions¶
A comprehensive presentation of the Glauber model can be found in e.g. https://inspirehep.net/record/742696. There are two main inputs to the Glauber model: the distribution of nucleons inside the nuclei, and the inelastic nucleon-nucleon cross-section. In very simplified terms, if the nuclei colliding in heavy ion collisions are taken as bags of nucleons, the Glauber model describes the collision of two nuclei as the incoherent collision of the all the nucleons inside the nuclei. This is illustrated in the following figure.
[Ref: Hannu Holopainen, PhD thesis ]
Every dark circle represents a nucleon that underwent a nucleon-nucleon collisions. They are referred to as “participant nucleons”. Light circles are nucleons that did not collide with any nucleons from the other nuclei. Nucleon-nucleon collisions are called binary collisions.
The Glauber model can be used as initial condition for the energy density \(\epsilon(X)\) in the plane transverse to the collision axis (x-y plane). Energy can be deposited either at the position of the binary nucleon collisions, or at the position of the participant nucleons, or a combination of the two. An input parameter in MUSIC controls binary_collision_scaling_fraction the fraction of energy density initialized with binary collision and participants nucleons.
The picture above is technically a Monte Carlo implementation of Glauber, referred to as Monte Carlo Glauber or MC Glauber. There is another version of Glauber, optical Glauber, which is roughly equivalent to averaging a large number of MC Glauber initial conditions, producing a smooth energy density distribution. Both MC Glauber and optical Glauber are available in MUSIC.
Here is a summary of which parameters are used for optical Glauber, which is used for MC Glauber, and which is used for both:
###################################
# parameters for initial conditions
###################################
Initial_profile 3 # type of initial condition
# 0: for Gubser flow test (future)
# 1: Optical Glauber model
# 3: MC-Glauber model
# 8: Read in initial profile from a file
# (e.g., IP-Glasma)
#
initialize_with_entropy 0 # 0: with energy density
# 1: with entropy density
#
s_factor 30.0 # normalization factor for the energy/entropy
# applied to MC-Glauber initial conditions
# (Initial_profile 3)
# or initial data file (Initial_profile 8)
#
Eta_plateau_size 20. # size of the plateau in spatial rapidity eta_s
Eta_fall_off 0.7 # the scale of the fall off of the
# plateau in eta_s direction
#
# parameters for Glauber model (Optical and MC Glauber, Initial_profile == 1 or 3)
binary_collision_scaling_fraction 0.14 # for wounded nucleon/binary
# collision mixing ratio in the
# Glauber model
SigmaNN 42.1 # NN inelastic cross section for
# Glauber model
Target Au # type of target nucleus
Projectile Au # type of projectile nucleus
# parameters for optical Glauber model only (effective only when Initial_profile == 1)
Impact_parameter 9. # impact parameter (fm)
#
Maximum_energy_density 54. # maximum energy density for Glauber
# initial condition
# (used instead of parameter "s_factor")
# parameters for MC Glauber model only (effective only when Initial_profile == 3)
sigma_0 0.4 # the width of the gaussian
# for entropy density profile
bmin 0.0 # minimum value of impact parameter
bmax 2.0 # maximum value of impact parameter
Note on the Glauber implementation in MUSIC
The implementation of MC Glauber used in MUSIC is described in Section IV of https://arxiv.org/pdf/1109.6289.pdf. This implementation may differ from other implementations of MC-Glauber, in particular in the following aspects:
- There is no minimum distance between nucleon sampled from the nucleus
- Negative binomial fluctuations are not included in the code, that is all nucleon-nucleon collisions deposit the same amount of energy/entropy
- Energy/entropy is deposited at the center of mass of binary collisions.
Since the Glauber model only provides an initialization for the energy density \(\epsilon(X)\) in the transverse plane, it must be complemented by other models or assumptions to provide a complete initial conditions. To obtain a profile in the \(\eta_s\), the input parameters Eta_plateau_size and Eta_fall_off can be used to multiply the transverse energy profile by a function that has a plateau around midrapidity (\(\eta_s \lesssim\) Eta_plateau_size ) and falls off with a width Eta_fall_off after.
This still only provides an initialization for the energy density \(\epsilon(X)\). When Glauber initial conditions are used, MUSIC initialises the flow velocity \(u^\mu(X)\) to \((1,0,0,0)\) and the viscous parts of the energy-momentum tensor (\(\pi^{\mu\nu}(X)\) and \(\Pi(X)\)) to zero. This can only be modified by modifying the code in “init.cpp”.