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\]

MUSIC can also solve for the conservation of the baryon density:

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

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 ‘Cell_small’ class:

class Cell_small {
 public:
        double epsilon = 0;
        double rhob    = 0;
        FlowVec u;

        ViscousVec Wmunu;
        double pi_b    = 0.;
};

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.

_images/contour_w_flow.png

[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 “eos*.cpp” files. 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. A range of equation of states are available within MUSIC .. In general, there shouldn’t be any need .. to modify the “eos*.cpp” files, since the equation of state which 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

    # EOS_to_use:
# 0: ideal gas
# 1: EOS-Q from azhydro
# 2: lattice EOS from Huovinen and Petreczky
# 3: lattice EOS from Huovinen and Petreczky
#    with partial chemical equilibrium (PCE) at 150 MeV
#    (see https://wiki.bnl.gov/TECHQM/index.php/QCD_Equation_of_State)
# 4: PCE EOS with chemical freeze out at 155 MeV
# 5: PCE EOS at 160 MeV
# 6: PCE EOS at 165 MeV
# 7: lattice EOS with CE match with UrQMD's hadron content
    # 8: lattice EOS from fits to the Wuppertal-Budapest lattice+HRG results
    # 9: lattice EOS obtained by matching HotQCD lattice results to the SMASH HRG ( see https://github.com/j-f-paquet/eos_maker  )
# 10-14: finite muB EOS from A. Monnai
    # 17: EOS from the BEST Collaboration

Warning

Some of the equation of states implemented in MUSIC are dated. At zero baryon chemical potentials,

EOS_to_use 9

should be a reasonable modern choice.

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.

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 “transport_coeffs.cpp” for the first order transport coefficients. Here is an abridged overview of these functions:

double TransportCoeffs::get_eta_over_s(double T) const {
        double eta_over_s;
        if (DATA.T_dependent_shear_to_s == 1) {
                eta_over_s = get_temperature_dependent_eta_over_s_default(T);
        } else if (DATA.T_dependent_bulk_to_s == 2) {
                eta_over_s = get_temperature_dependent_eta_over_s_duke(T);
        } else if (DATA.T_dependent_bulk_to_s == 3) {
                eta_over_s = get_temperature_dependent_eta_over_s_sims(T);
        } else {
                eta_over_s = DATA.shear_to_s;
        }
        return eta_over_s;
}


double TransportCoeffs::get_temperature_dependent_eta_over_s_default(double T) const {
        double Ttr = 0.18/hbarc;  // phase transition temperature
        double Tfrac = T/Ttr;
        double eta_over_s;
        if (T < Ttr) {
                eta_over_s = (DATA.shear_to_s + 0.0594*(1. - Tfrac)
                                          + 0.544*(1. - Tfrac*Tfrac));
        } else {
                eta_over_s = (DATA.shear_to_s + 0.288*(Tfrac - 1.)
                                          + 0.0818*(Tfrac*Tfrac - 1.));
        }
        return(eta_over_s);
}


double TransportCoeffs::get_temperature_dependent_eta_over_s_duke(
                                                                                                double T_in_fm) const {
        double T_in_GeV = T_in_fm*hbarc;
        double Ttr_in_GeV = 0.154;
        double Tfrac = T_in_GeV/Ttr_in_GeV;

        double eta_over_s = (DATA.eta_over_s_min
                                                 + (DATA.eta_over_s_slope)*(T_in_GeV - Ttr_in_GeV)
                                                   *pow(Tfrac,DATA.eta_over_s_curv));
        return eta_over_s;
}


double TransportCoeffs::get_temperature_dependent_eta_over_s_sims(
                                                                                                double T_in_fm) const {
        double T_in_GeV = T_in_fm*hbarc;
        double T_kink_in_GeV = DATA.eta_over_s_T_kink_in_GeV;
        double low_T_slope = DATA.eta_over_s_low_T_slope_in_GeV;
        double high_T_slope = DATA.eta_over_s_high_T_slope_in_GeV;
        double eta_over_s_at_kink = DATA.eta_over_s_at_kink;

        const double eta_over_s_min = 1e-6;
        double eta_over_s;
        if (T_in_GeV < T_kink_in_GeV) {
                eta_over_s = (eta_over_s_at_kink
                                          + low_T_slope*(T_in_GeV - T_kink_in_GeV));
        } else {
                eta_over_s = (eta_over_s_at_kink
                                          + high_T_slope*(T_in_GeV - T_kink_in_GeV));
        }
        eta_over_s = std::max(eta_over_s,eta_over_s_min);
        return(eta_over_s);
}

The normalization of the shear and bulk relaxation times can also be adjusted with input parameters.

TransportCoeffs::TransportCoeffs(const EOS &eosIn, const InitData &Data_in)
        : DATA(Data_in), eos(eosIn) {
        shear_relax_time_factor_ = DATA.shear_relax_time_factor;
        bulk_relax_time_factor_  = DATA.bulk_relax_time_factor;
}

The other second-order transport coefficients are set in “transport_coeffs.h”:

// shear
double get_tau_pipi_coeff() const { return(10./7.); }
double get_delta_pipi_coeff() const { return(4./3.); }
double get_phi7_coeff() const { return(9./70.); }
double get_lambda_piPi_coeff() const { return(6./5.); }

// bulk
double get_lambda_Pipi_coeff() const { return(8./5.); }
double get_delta_PiPi_coeff() const { return(2./3.); }
double get_tau_PiPi_coeff() const { return(0.); }  // unknown

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. In fact, in heavy ion collisions, particlization based on Cooper-Frye requires that hadronization occurs 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.

_images/contour_w_flow_particlization.png

[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 files 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” files are then used the evaluate “Cooper-Frye” when MUSIC is run in mode 3.