2. Hydrodynamics

Hydrodynamics is used to describe the space-time evolution of the energy-momentum tensor \(T^{\mu\nu}\) of the strongly-coupled quark-gluon plasma.

This makes possible to learn about certain general many-body properties of QCD without necessarily tracking the complicated dynamics of every quark, gluon and hadron in the system. Such properties of QCD include the equation of state and the shear and bulk viscosities of QCD.

The basis of hydrodynamics is local conservation of energy, which can be written:

\[\partial_\mu T^{\mu\nu}(X)=0\]

Additional conserved charges

If there are conserved charges in the fluid, these are written

\[\partial_\mu j^\mu(X) =0\]

The public version of MUSIC does not have any charge conservation and as such, this equation is not solved in the code.

The energy-momentum tensor \(T^{\mu\nu}\) can be decomposed into four different fields:

\[T^{\mu \nu}(X)=\epsilon(X) u^\mu(X) u^\nu(X)- \left[ \mathcal{P}(X)+\Pi(X) \right] \left[ g^{\mu\nu}-u^\mu(X) u^\nu(X) \right ] +\pi^{\mu\nu}(X)\]

where

• \(\epsilon(X)\) is the energy density
• \(u^\mu(X)\) is the flow velocity
• \(\pi^{\mu\nu}(X)\) is the shear stress tensor
• \(\Pi(X)\) is the bulk pressure

The fifth field, \(\mathcal{P}(X)\) , is the pressure and it is not independent: it is related to the energy density by the equation of state, as we will discuss below.

The tensor \(g^{\mu\nu}\) is the metric tensor.

Warning

Beware of the metric convention. These notes are based on the mostly-minus convention:

\(g^{\mu\nu}=diag(1,-1,-1,-1)\)

The metric convention in most if not all MUSIC publications is mostly-minus.

On the other hand, the formulae implemented in the code use the mostly-plus metric. Unless you have to modify the core algorithms of MUSIC, this should not affect you in any way. If you do need to modify the code, pay special attention to the metric.

What the decomposition of \(T^{\mu\nu}\) above means is that every spacetime point of the quark-gluon plasma has a value for \(T^{\mu\nu}\), which can be expressed in terms of \(\epsilon\), \(u^\mu\), \(\pi^{\mu\nu}\) and \(\Pi\). In MUSIC, this information is saved in the ‘Grid’ class:

class Grid
{
 public:
    double epsilon;
    double p;
    [...]

    //! u[flag][mu]: flag=0 is the actual values. flag !=0 are for RK steps
    double **u;

    [...]
    //! bulk part of the TJb with the rk_flag
    double ***Pimunu;
    [...]

    //! bulk pressure */
    double *pi_b;

    [...]
};

Note that the “Grid” class contains much more information, which was hidden here and replaced by “[…]”. Most of this information is used to store intermediate numerical steps. Other is redundant information, for example the temperature corresponding to the energy density, as calculated with the equation of state.

2.1. Ideal hydrodynamics

In ideal hydrodynamics, \(\pi^{\mu\nu}\) and \(\Pi\) are zero, and the energy-momentum tensor reduces to:

\[T^{\mu \nu}(X)=\epsilon(X) u^\mu(X) u^\nu(X)- \mathcal{P}(X) \left[ g^{\mu\nu}-u^\mu(X) u^\nu(X) \right ]\]

The flow velocity \(u^\mu(X)\) is defined as the unit vector satisfying:

\[T^\mu_{\phantom{\mu}\nu}(X) u^\nu(X) \equiv \epsilon(X) u^\mu(X); \,\,u^\mu(X) u_\mu(X)=1\]

meaning that the flow velocity represents the flow of energy. This is the so-called “Landau frame”. Other choices of frames are possible, but the Landau frame is the one used in MUSIC.

Because it is a normalised four-vector, \(u^\mu(X)\) has three degrees of freedom for each spacetime point. The energy density \(\epsilon(X)\) is another degree of freedom.

Since the pressure \(\mathcal{P}(X)\) is related to the energy density, \(T^{\mu\nu}(X)\) has four degrees of freedom per spacetime point, and the energy conservation equation \(\partial_\mu T^{\mu\nu}(X)=0\) provides four constraints per spacetime point: this forms a closed set of differential equations.

[Ref.: J.-F.P.]

MUSIC can be run ideal hydrodynamics by turning off both shear and bulk viscosities in the input file:

Include_Shear_Visc_Yes_1_No_0 0    # include shear viscous effect
Include_Bulk_Visc_Yes_1_No_0 0     # include bulk viscous effect
Include_second_order_terms 0       # include second order non-linear coupling terms

Pressure, equation of state and the assumptions behind ideal hydrodynamics

The pressure is related to the energy density because ideal hydrodynamics is built around the assumption of local thermal equilibrium: one can define a small region at each spacetime point X where the fluid is assumed to be in local equilibrium. The system evolves in space and time, but the fluid is assumed to always maintain equilibrium locally.

Thus one can use thermodynamics locally to relate the energy density to the pressure, the temperature, the entropy: \(\mathcal{P}(X)=\mathcal{P}(\epsilon(X))\) , …

The equation of state is that of the strongly-coupled quark-gluon plasma. It is calculated with lattice QCD and matched to a hadron resonance gas [*].

[*]	The composition of this hadron resonance gas is important at the “particlization” stage, when the energy-momentum tensor of the fluid is converted into hadrons. This discussion is best postponed to the “Particlization” section.

Since the equation of state calculated from lattice QCD and from the hadron resonance gas are numerical (no simple closed expression), and since the matching between the two may also be performed numerically, the equation of state must effectively be provided to MUSIC either in the form of a parametrization or a table. Tables are used for all current equations of state in MUSIC.

The code related to the equations of state is in the file “eos.cpp”. There are multiple functions to allocate and load to memory the files containing the equation of state information. There are also functions to interpolate the equation of state tables. In general, there shouldn’t be any need to modify “eos.cpp”, since the equation of state can be chosen from the input file with the parameter EOS_to_use.

The different choices of equations of state are explained at various points in the code and in example input files. Here is a summary:

EOS_to_use 2                  # type of the equation of state
                              # 0: ideal gas
                              # 1: EOS-Q from azhydro
                              # 2: lattice EOS s95p-v1
                              #    (from Huovinen and Petreczky)
                              # 3: lattice EOS s95p with partial
                              #    chemical equilibrium (PCE) at 150 MeV
                              #    (see https://wiki.bnl.gov/TECHQM
                              #         /index.php/QCD_Equation_of_State)
                              # 4: lattice EOS s95p with chemical freeze
                              #    out at 155 MeV
                              # 5: lattice EOS s95p at 160 MeV
                              # 6: lattice EOS s95p at 165 MeV
                              # 7: lattice EOS s95p-v1.2 for UrQMD

Warning

Many of the equation of states implemented in MUSIC are dated. Updated ones are being implemented and tested: they should be available soon.

2.2. Viscous hydrodynamics

Almost all hydrodynamics models of heavy ion collisions use Israel-Stewart-type second-order viscous hydrodynamics.

In these models, conservation of energy \(\partial_\mu T^{\mu\nu}(X)=0\) and the equation of state \(\mathcal{P}=\mathcal{P}(\epsilon)\) are complemented by equations of motion for the shear stress tensor and the bulk pressure.

The equation of motions implemented in the public version of MUSIC are based on these publications:

• Denicol, G. S.; Niemi, H.; Molnár, E.; Rischke, D. H. “Derivation of transient relativistic fluid dynamics from the Boltzmann equation”, PRD85:114047 (2012); PRD91:039902 (2015)
• Molnár, E.; Niemi, H.; Denicol, G. S.; Rischke, D. H. “Relative importance of second-order terms in relativistic dissipative fluid dynamics” PRD89:074010 (2014)

This derivation of second-order viscous hydrodynamics is sometimes referred to as DNMR, after its authors.

The equations of motion are also summarised in this publication [†]:

• Denicol, G. S.; Jeon, S.; Gale, C. “Transport coefficients of bulk viscous pressure in the 14-moment approximation” PRC90:024912 (2014)

[†]

There is a \(\tau_{\pi}\) missing in front of \(\pi _{\alpha }^{\left\langle \mu \right.}\omega ^{\left. \nu \right\rangle \alpha }\) in Equation 2 of “Denicol, G. S.; Jeon, S.; Gale, C. PRC90:024912 (2014)” This is clarified in the erratum of “Denicol, G. S.; Niemi, H.; Molnár, E.; Rischke, D. H. PRD85:114047 (2012)” : PRD91:039902 (2015)

Explicitly, the equations of motions can be compactly written as:

\[\begin{split}\tau _{\pi }\dot{\pi}^{\left\langle \mu \nu \right\rangle }+\pi ^{\mu \nu } &=&2\eta \sigma ^{\mu \nu }+2 \tau_{\pi} \pi _{\alpha }^{\left\langle \mu \right. }\omega ^{\left. \nu \right\rangle \alpha }-\delta _{\pi \pi }\pi ^{\mu \nu }\theta +\varphi _{7}\pi _{\alpha }^{\left\langle \mu \right. }\pi ^{\left. \nu \right\rangle \alpha }-\tau _{\pi \pi }\pi _{\alpha }^{\left\langle \mu \right. }\sigma ^{\left. \nu \right\rangle \alpha } \notag \\ &&\;+\lambda _{\pi \Pi }\Pi \sigma ^{\mu \nu }+\varphi _{6}\Pi \pi ^{\mu \nu}\end{split}\]

for the shear stress tensor and

\[\tau _{\Pi }\dot{\Pi}+\Pi =-\zeta \theta -\delta _{\Pi \Pi }\Pi \theta +\varphi _{1}\Pi ^{2}+\lambda _{\Pi \pi }\pi ^{\mu \nu }\sigma _{\mu \nu }+\varphi _{3}\pi ^{\mu \nu }\pi _{\mu \nu }\]

for the bulk pressure. Please refer to “Denicol, G. S.; Jeon, S.; Gale, C. PRC90:024912 (2014)” for more information about the notation used to write the equation of motions.

Vorticity

Although the vorticity term, \(2 \tau_{\pi} \pi _{\alpha }^{\left\langle \mu \right.}\omega ^{\left. \nu \right\rangle \alpha }\) is included in the code, it is turned off by default. Turning it on requires modifying the following variable in the code to “1” instead of “0”:

int include_Vorticity_term = 0;

in function “Diss::Make_uWSource()” of “dissipative.cpp”

These equations have two first-order transport coefficients, the shear viscosity \(\eta\) and the bulk viscosity \(\zeta\).

There are also 11 second-order transport coefficients:

\[ \begin{align}\begin{aligned}\tau _{\pi }, \delta _{\pi \pi }, \varphi _{7}, \tau _{\pi \pi }, \lambda _{\pi \Pi }, \varphi _{6}\\\tau _{\Pi }, \delta _{\Pi \Pi }, \varphi _{1},\lambda _{\Pi \pi }, \varphi _{3}\end{aligned}\end{align} \]

All these transport coefficients, first and second order, are characteristic properties of how the quark-gluon plasma responds to deviation from equilibrium.

Transport coefficients from first principles

In principle, the transport coefficients of the quark-gluon plasma are calculable from QCD. In practice, these are extremely challenging calculations which can only be performed under very specific conditions, such as very high temperature (possibly temperatures orders of magnitude higher than the typical temperatures of the quark-gluon plasma) or with effective models of QCD.

In MUSIC, the second order transport coefficients are related to the first-order transport coefficients through relations derived in the three publications listed above:

\[\begin{split}\tau_\pi & = & \frac{5 \eta}{(\epsilon+\mathcal{P})} \\ \delta_{\pi\pi} & = & \frac{4}{3} \tau_\pi \\ \phi_7 & = & \frac{9}{70 \mathcal{P}} \\ \tau_{\pi\pi} & = & \frac{10}{7} \tau_\pi \\ \lambda_{\pi\Pi} & = & \frac{6}{5} \\ \tau_{\Pi} &=& \frac{\zeta}{15 \left( \frac{1}{3}-c_s^2 \right)^2\left( \epsilon+\mathcal{P} \right)} \\ \delta_{\Pi\Pi} &=& \frac{2}{3} \tau_{\Pi} \\ \lambda_{\Pi\pi} &=& \frac{8}{5} \left( \frac{1}{3}-c_s^2 \right) \tau_{\Pi} \\\end{split}\]

Three second order transport coefficient have not been derived in the publications above and are set to zero: \(\varphi _{1}=\varphi _{3}=\varphi _{6}=0\).

The first order transport coefficients are parametrized, with the objectives that they will be extracted from measurements.

The relevant part of the code can be found in functions “Diss::Make_uWSource()” and “Diss::Make_uPiSource()” in source file “dissipative.cpp”. Here is an abridged overview of these functions:

double Diss::Make_uWSource(...) {

    [...]

    if (DATA->T_dependent_shear_to_s == 1) {
        if (T < Ttr) {
            shear_to_s=0.681-0.0594*T/Ttr-0.544*(T/Ttr)*(T/Ttr);
        } else {
            shear_to_s=-0.289+0.288*T/Ttr+0.0818*(T/Ttr)*(T/Ttr);
        }
    } else {
        shear_to_s = DATA->shear_to_s;
    }

    int include_WWterm = 0;
    int include_Vorticity_term = 0;
    int include_Wsigma_term = 0;
    if (DATA->include_second_order_terms == 1) {
        include_WWterm = 1;
        include_Wsigma_term = 1;
    }
    if (DATA->Initial_profile == 0) {
        include_WWterm = 0;
        include_Wsigma_term = 0;
        include_Vorticity_term = 0;
    }

// ////////////////////////////////////////////////////////////////////// //
// ////////////////////////////////////////////////////////////////////// //
//                 Defining transport coefficients                        //
// ////////////////////////////////////////////////////////////////////// //
// ////////////////////////////////////////////////////////////////////// //
 //double s_den = eos->get_entropy(epsilon, rhob);
 shear = (shear_to_s)*(grid_pt->epsilon + grid_pt->p)/(T + 1e-15);
 tau_pi = 5.0*shear/(grid_pt->epsilon + grid_pt->p + 1e-15);

[...]

 // transport coefficient for nonlinear terms -- shear only terms -- 4Mar2013
 double transport_coefficient, transport_coefficient2, transport_coefficient3, transport_coefficient4;
 // transport coefficients of a massless gas of single component particles
 transport_coefficient  = 9./70.*tau_pi/shear*(4./5.) ;
 transport_coefficient2 = 4./3.*tau_pi;
 transport_coefficient3 = 10./7.*tau_pi;
 transport_coefficient4 = 2.*tau_pi;

 // transport coefficient for nonlinear terms -- coupling to bulk viscous pressure -- 4Mar2013
 double transport_coefficient_b, transport_coefficient2_b;
 // transport coefficients not yet known -- fixed to zero
 transport_coefficient_b  = 6./5.*tau_pi ;
 transport_coefficient2_b = 0.;


/* This source has many terms */
/* everything in the 1/(tau_pi) piece is here */
/* third step in the split-operator time evol
   use Wmunu[rk_flag] and u[rk_flag] with rk_flag = 0 */

// ////////////////////////////////////////////////////////////////////// ///
// ////////////////////////////////////////////////////////////////////// ///
//            Wmunu + transport_coefficient2*Wmunu*theta                  ///
// ////////////////////////////////////////////////////////////////////// ///
// ////////////////////////////////////////////////////////////////////// ///

// full term is
    tempf = -(1.0 + transport_coefficient2*(grid_pt->theta_u[rk_flag]) )*(grid_pt->Wmunu[rk_flag][mu][nu]);
// FOR GUBSER ANALYTIC
//     tempf = -( transport_coefficient2*(grid_pt->theta_u[rk_flag]) )*(grid_pt->Wmunu[rk_flag][mu][nu]);


// ////////////////////////////////////////////////////////////////////// ///
// ////////////////////////////////////////////////////////////////////// ///
//             Navier-Stokes Term -- -2.*shear*sigma^munu                 ///
// ////////////////////////////////////////////////////////////////////// ///
// ////////////////////////////////////////////////////////////////////// ///

// remember: dUsup[m][n] = partial^n u^m  ///
// remember:  a[n]  =  u^m*partial_m u^n  ///

    for( a=0;a<4;a++ )
    {
        for( b=0;b<4;b++ )
        {
            sigma[a][b] = grid_pt->sigma[rk_flag][a][b];
        }
    }

[...]

}/* Make_uWSource */


double Diss::Make_uPiSource([...]) {

    [...]

    // switch to include non-linear coupling terms in the bulk pi evolution
    int include_BBterm = 0;
    int include_coupling_to_shear = 0;
    if (DATA->include_second_order_terms == 1) {
        include_BBterm = 1;
        include_coupling_to_shear = 1;
    }

    [...]

    // T dependent bulk viscosity from Gabriel
    // ///////////////////////////////////////////
    //            Parametrization 1             //
    // ///////////////////////////////////////////
    double Ttr=0.18/0.1973;
    double dummy=temperature/Ttr;
    double A1=-13.77, A2=27.55, A3=13.45;
    double lambda1=0.9, lambda2=0.25, lambda3=0.9, lambda4=0.22;
    double sigma1=0.025, sigma2=0.13, sigma3=0.0025, sigma4=0.022;

    bulk = A1*dummy*dummy + A2*dummy - A3;
    if(temperature < 0.995*Ttr)
    {
        bulk = lambda3*exp((dummy-1)/sigma3)+ lambda4*exp((dummy-1)/sigma4)+0.03;
    }
    if(temperature > 1.05*Ttr)
    {
        bulk = lambda1*exp(-(dummy-1)/sigma1)+ lambda2*exp(-(dummy-1)/sigma2)+0.001;
    }
    bulk = bulk*(grid_pt->epsilon + grid_pt->p)/temperature;

    [...]

    // defining bulk relaxation time and additional transport coefficients
    // Bulk relaxation time from kinetic theory
    Bulk_Relax_time    = 1./14.55/(1./3.-cs2)/(1./3.-cs2)/(grid_pt->epsilon + grid_pt->p)*bulk;

    transport_coeff1   = 2.0/3.0*(Bulk_Relax_time);          // from kinetic theory, small mass limit
    transport_coeff2   = 0.;                                 // not known; put 0
    transport_coeff1_s = 8./5.*(1./3.-cs2)*Bulk_Relax_time;  // from kinetic theory
    transport_coeff2_s = 0.;                                 // not known;  put 0

    // Computing Navier-Stokes term (-bulk viscosity * theta)
    NS_term = -bulk*(grid_pt->theta_u[rk_flag]);

    // Computing relaxation term and nonlinear term: - Bulk - transport_coeff1*Bulk*theta
    tempf = -(grid_pt->pi_b[rk_flag]) - transport_coeff1*(grid_pt->theta_u[rk_flag])*(grid_pt->pi_b[rk_flag]);

    // Computing nonlinear term: + transport_coeff2*Bulk*Bulk
    if (include_BBterm == 1)
        BB_term = transport_coeff2*(grid_pt->pi_b[rk_flag])*(grid_pt->pi_b[rk_flag]);
    else
        BB_term = 0.0;

    // Computing terms that Couple with shear-stress tensor
    double Wsigma, WW, Shear_Sigma_term, Shear_Shear_term, Coupling_to_Shear;

    if (include_coupling_to_shear == 1)
    {
        // Computing sigma^mu^nu
        double sigma[4][4];
        for(int a=0;a<4;a++ )
        {
            for(int b=0;b<4;b++ )
            {
                sigma[a][b] = grid_pt->sigma[rk_flag][a][b];
            }
        }

        Wsigma = (
                    grid_pt->Wmunu[rk_flag][0][0]*sigma[0][0]
                   +grid_pt->Wmunu[rk_flag][1][1]*sigma[1][1]
                   +grid_pt->Wmunu[rk_flag][2][2]*sigma[2][2]
                   +grid_pt->Wmunu[rk_flag][3][3]*sigma[3][3]

                   -2.*(
                         grid_pt->Wmunu[rk_flag][0][1]*sigma[0][1]
                        +grid_pt->Wmunu[rk_flag][0][2]*sigma[0][2]
                        +grid_pt->Wmunu[rk_flag][0][3]*sigma[0][3]
                        )

                   +2.*(
                         grid_pt->Wmunu[rk_flag][1][2]*sigma[1][2]
                        +grid_pt->Wmunu[rk_flag][1][3]*sigma[1][3]
                        +grid_pt->Wmunu[rk_flag][2][3]*sigma[2][3]
                        )
                 );

        [...]

        // multiply term by its respective transport coefficient
        Shear_Sigma_term = Wsigma*transport_coeff1_s;

        [...]

    }
    else
        Coupling_to_Shear = 0.0;

    [...]

}/* Make_uPiSource */

Note that the bulk viscosity is hard-coded. For the shear viscosity, there are two options: a hard-coded parametrization or a constant value of \(\eta/s\), whose value is set by the input parameter Shear_to_S_ratio. The hard-coded parametrizations for \(\eta/s\) and \(\zeta/s\) are shown below:

The input parameters that control shear and bulk viscosity are:

# transport coefficients
Viscosity_Flag_Yes_1_No_0 1        # turn on viscosity in the evolution

Include_Shear_Visc_Yes_1_No_0 1    # include shear viscous effect
Shear_to_S_ratio 0.08              # value of \eta/s
T_dependent_Shear_to_S_ratio  0    # flag to use temperature dep. \eta/s(T)

Include_Bulk_Visc_Yes_1_No_0 0     # include bulk viscous effect
Include_second_order_terms 0       # include second order non-linear coupling terms

The Include_second_order_terms input parameter turns on and off a subset of second-order terms: \(\varphi _{7}\pi _{\alpha }^{\left\langle \mu \right. }\pi ^{\left. \nu \right\rangle \alpha }\), \(\tau _{\pi \pi }\pi _{\alpha }^{\left\langle \mu \right. }\sigma ^{\left. \nu \right\rangle \alpha }\), \(\lambda _{\pi \Pi }\Pi \sigma ^{\mu \nu }\) and \(\varphi _{6}\Pi \pi ^{\mu \nu}\) in the shear tensor equation of motion, and \(\delta _{\Pi \Pi }\Pi \theta\), \(\lambda _{\Pi \pi }\pi ^{\mu \nu }\sigma _{\mu \nu}\) and \(\varphi _{3}\pi ^{\mu \nu }\pi _{\mu \nu }\) in the bulk pressure equation of motion.

2.3. Solving the hydrodynamic equations: Numerical parameters

The previous section was a brief overview of the physics and the equations that are solved in MUSIC. Ultimately the equations of hydrodynamics are a set of coupled differential equations that must be solved numerically.

Hyperbolic coordinates are used in MUSIC:

\[\left(\tau=\sqrt{t^2-z^2},x,y,\eta_s=\textrm{arctanh}(z/t)\right)\]

The “spatial” coordinates \((x,y,\eta_s)\) are discretized on a grid, and the “time” coordinate \(\tau\) is evolved. Having a grid in \((x,y,\eta_s)\) implies that each grid points have a coordinates which must be specified in one way or another. The grid used in MUSIC is regular in \((x,y,\eta_s)\). It is centered around the origin, \((x,y,\eta_s)=\vec{0}\), up to a small detail.

Coordinate system of the grid used in MUSIC

The \((x,y,\eta_s)\) grid in MUSIC is centered at \(x=0\) and \(y=0\) in the sense that there is an equal number of grid points on each side of \(x=0\) and \(y=0\). This is achieved by having an odd number of grid points in x and y. The number of grid points specified by the user is even, but an extra grid point is added so that the total number of grid points is odd and there can be a grid point exactly at \(x=y=0\).

The number of grid point in \(\eta_s\) specified by the user is even, and no additional grid point is added by MUSIC. To make sure that there is a grid point at \(\eta_s=0\), the grid is staggered to the left: there is one more grid point in \(\eta_s<0\) than there is for \(\eta_s>0\).

There only time these subtleties should be of concern is when external initial conditions are loaded into MUSIC.

Since it is a regular grid, the grid’s coordinates can be specified by the grid’s total size in each direction (\(x/y/\eta_s\)) and the number of cells that each coordinate is divided in.

For the \(\tau\) coordinate, an initial time is specified, along with a maximum time and a time-step.

The corresponding input parameters are (see also this file):

Initial_time_tau_0 0.4        # starting time of the hydrodynamic
                              # evolution (fm/c)

Delta_Tau 0.02                # time step to use in the evolution [fm/c]

Total_evolution_time_tau 300.  # the maximum allowed running
                              # evolution time (fm/c)
                              # almost always needs to be set to some large number

Eta_grid_size 12.8            # spatial rapidity range
                              # [-Eta_grid_size/2, Eta_grid_size/2 - delta_eta]
Grid_size_in_eta 64           # number of the grid points in spatial
                              # rapidity direction
                              # Must have at least 4 cells per processor.
                              # Must be an even number.
                              # One cell is positioned at eta=0,
                              # half the cells are at negative eta,

                              # the rest (one fewer) are at positive eta
X_grid_size_in_fm 20.0        # spatial range along x direction in the
                              # transverse plane
                              # [-X_grid_size_in_fm/2, X_grid_size_in_fm/2]
Y_grid_size_in_fm 20.0        # spatial range along y direction in the
                              # transverse plane
                              # [-Y_grid_size_in_fm/2, Y_grid_size_in_fm/2]
Grid_size_in_y 200            # number of the grid points in y direction
Grid_size_in_x 200            # number of the grid points in x direction

The algorithm used in MUSIC to solve the equations of hydrodynamics is named Kurganov-Tadmor. It is an algorithm designed to handle large gradients when solving partial differential equations. Kurganov-Tadmor is itself a variation of a larger family of partial differential equation solver called “MUSCL” (MUSIC stands for “MUSCL-based Ion Collider”).

A small number of parameters control numerical details of how the equations are solved. These parameters are hard-coded in the public version of MUSIC, and you should not need to modify them unless you are an expert with numerical solution of hydrodynamics with MUSCL-based schemes.

2.4. Construction of the particlization hypersurface

Hydrodynamics cannot describe a fluid of arbitrary low density, in particular if this fluid is expanding rapidly (as is the case in heavy ion collisions). Physically, this means that the hydrodynamics description of the plasma produced in heavy ion collisions must be stopped at some point, and that the matter being described with hydrodynamics must be converted into different degrees of freedom. This procedure has been given different names over the years, but “particlization” is probably the most accurate and descriptive, since the fluid is converted into particles. This conversion from fluid to particle is performed in MUSIC using the Cooper-Frye algorithm, as is the case for almost every simulation of heavy ion collisions.

Particlization, Cooper-Frye, hadronization and freeze-out

“Particlization” is not necessarily the most widely used term to describe the conversion of the fluid to particles. “Cooper-Frye” is probably used more often.

It is important not to confuse “particlization/Cooper-Frye” with “hadronization” and “freeze-out”, which are sometimes used loosely to refer to particlization.

Technically speaking, “hadronization” refers to the reconfinement of the strongly-coupled quark-gluon plasma into hadrons. However, hadronization and particlization do not have to happen at the same time. Cooper-Frye actually requires that hadronization happens before particlization.

“Freeze-out” refers to the moment when interactions among particles stop completely (or, in practice, when these interactions become negligible). For a hadron gas, it implies that the gas is so dilute that hadrons do not interact with each other anymore. “Freeze-out” can also refer to two separate concepts: chemical freeze-out, when inelastic collisions stop, and kinetic freeze-out, when both inelastic and elastic collisions stop. When “freeze-out” is used alone, it usually refers to the kinetic freeze-out.

The Cooper-Frye formula will be discussed in a later section. What is relevant to discuss here is the particlization criteria, the criteria that is used to stop the hydrodynamics and convert the fluid into hadrons. In MUSIC as in almost all other hydrodynamics simulation of heavy ion collisions, particlization is assumed to happen at a given temperature (or energy density, which is related to temperature through the equation of state). As the temperature of the hydrodynamic medium decreases, different parts of the medium reach the particlization temperature. Connecting all the points at which the local temperature of the medium cross the particlization temperature forms a 4D spacetime hypersurface. This hypersurface is complete once every point in the hydrodynamic medium is below this particlization criteria, and hydrodynamics can then be stopped.

[Ref.: J.-F.P.]

A two-dimensional view of a particlization hypersurface is shown on the figure above. The colour and the contours represent different temperatures. The horizontal axis is time and the vertical axis is in the transverse direction of the fluid (perpendicular to the collisions axis). The particlization criteria is \(T=150\) MeV in this example.

Warning

All particlization parameters in MUSIC are labelled with “freeze” or “freeze_out”, since the initial studies performed with MUSIC did not distinguish between particlization and freeze-out.

In MUSIC, this hypersurface is constructed at every timestep in \(\tau\) of the numerical solver. There is a number of parameters controlling this feature of MUSIC:

Do_FreezeOut_Yes_1_No_0 1       # Whether the particlization hypersurface is calculated or not
                                # If the hypersurface is not calculated, the code will evolve
                                # in time :math:`\tau` until the maximum time specified by
                                # input parameter "Total_evolution_time_tau"

average_surface_over_this_many_time_steps 5   # Coarse-grain the hypersurface in the tau

use_eps_for_freeze_out 1        # flag to use energy density as criteria to
                                # find particlization surface
                                # 0: use temperature, 1: use energy density
T_freeze 0.135                  # particlization temperature (GeV)
epsilon_freeze 0.18             # the particlization energy density (GeV/fm^3)

The hypersurface is saved in a file named “surface.dat”. The information that is saved is the vector normal to the hypersurface pointing outward from the fluid, the position of every hypersurface element, along with the value of all hydrodynamic fields at this spacetime point.

This “surface.dat” file is then used the evaluate “Cooper-Frye” when MUSIC is run in mode 3.