Optimized fluid dynamics for heavy ion collisions

TL;DR

This paper discusses the development and application of viscous anisotropic hydrodynamics for simulating the quark-gluon plasma in heavy ion collisions, showing improved performance over traditional viscous hydrodynamics in complex 3D scenarios.

Contribution

It introduces an efficient 3D simulation framework for viscous anisotropic hydrodynamics and compares its performance to accelerated viscous hydrodynamics simulations.

Findings

Anisotropic hydrodynamics outperforms viscous hydrodynamics in complex scenarios.

GPU-accelerated simulations enable faster 3D modeling.

Results demonstrate improved accuracy in describing quark-gluon plasma evolution.

Abstract

The (viscous) anisotropic hydrodynamic approach, especially after perturbative inclusion of all residual viscous terms, has been shown to dramatically outperform viscous hydrodynamics in several simplified situations for which exact solutions exist but which share with realistic expansion scenarios the problem of large dissipative currents. We will report on the present status of applying viscous anisotropic hydrodynamics in a highly efficient simulation of the full three-dimensional quark-gluon plasma. Results from accelerated $3{+}1$-dimensional viscous hydrodynamic simulations using graphics processing units will be compared to the anisotropic frameworks.

Figures

Tables

Table 1. Table 1: Performance results of the C/CPU and CUDA/GPU versions of CPU-VH and GPU-VH by measuring the computer time it takes to complete one full RK step, averaged over 100 time steps, at different spatial resolutions. v3.

Number of grid points	C/CPU	CUDA/GPU	Speedup
	(ms/step)	(ms/step)
128 $\times$ 128 $\times$ 32	7145.978	63.261	112.960
128 $\times$ 128 $\times$ 64	13937.896	123.527	112.833
128 $\times$ 128 $\times$ 128	30717.367	244.450	125.659
256 $\times$ 256 $\times$ 32	25934.547	236.593	109.617
256 $\times$ 256 $\times$ 64	57387.141	472.391	121.482
256 $\times$ 256 $\times$ 128	129239.959	939.340	137.586

Equations

$$T^{\mu \nu }=\mathcal{E}{u}^{\mu }{u}^{\nu }-{\mathcal{P}}_{\perp }{\Delta }_{\perp }^{\mu \nu }+{\hat{\mathcal{P}}}_L{z}^{\mu }{z}^{\nu }+2{\tilde{W}}_{\perp z}^{(\mu }{z}^{\nu )}+{\tilde{\pi }}_{\perp }^{\mu \nu }.$$

$${\hat{\mathcal{P}}}_{\perp }\equiv -\frac{1}{2}{\Delta }_{\perp }^{\mu \nu }{\hat{T}}_{\mu \nu }=\frac{1}{2}\left(\mathcal{E}-{\hat{\mathcal{P}}}_L-\bar{\mathcal{R}}(\xi )(\mathcal{E}-3{\mathcal{P}}_0)\right),$$

$$d\tilde{\Pi }$$

$$d{\hat{\mathcal{P}}}_L$$

$$d{\tilde{W}}_{\perp z}^{\mu }$$

$$d{\tilde{\pi }}_{\perp }^{\mu \nu }$$

$$u_{\perp }=\frac{{\beta }_{\perp }}{{\beta }_L+{\mathcal{P}}_{\perp }(u_{\perp })}\sqrt{1+u_{\perp }^2}.$$

$$u^0=\frac{1}{x}\sqrt{1+u_{\perp }^2},u^{1,2}=u_{\perp }\frac{{\bar{M}}^{1,2}}{{\bar{M}}_{\perp }},u^3=\frac{\tilde{F}}{\tau x}\sqrt{1+u_{\perp }^2}.$$

Full text

111Supported by the US Department of Energy under awards DE-SC0004286 and DE-SC0013470 and by the National Science Foundation through the JETSCAPE Collaboration (award number ACI-1550223).

1 Introduction

Quantitative modeling of relativistic heavy ion collisions has dramatically improved in recent years. Simulations include an initial state and pre-equilibrium model for the first ${\sim\,}0.5-1$ fm/$c$ when the quark-gluon plasma (QGP) is highly momentum anisotropic and far from thermalized, followed by viscous hydrodynamic simulations of the medium for the semi-isotropic and semi-thermal QGP for the next ${\cal O}(10)$ fm/$c$, and finally a microscopic hadronic treatment of the reformed hadrons. One of the main uncertainties in this chain is from the early evolution stage when the QGP is highly anisotropic and cannot be treated with viscous hydrodynamics. This uncertainty is usually encapsulated in model parameters describing the initial state and early pre-equilibrium dynamics. To constrain these and other model parameters with the help of experimental observables, these simulations are coupled with advanced statistical methods based on Bayesian statistics [1, 2]. Since it is very computationally expensive to train these models to cover very high-dimensional parameter spaces we optimize the simulations along two fronts: we (1) speed up the $3{+}1$ dimensional fluid dynamic simulation by performing the calculations on graphics processing units (GPUs) [3] and (2) limit the uncertainty from the early evolution stage by using the viscous anisotropic hydrodynamic formalism [4] which accounts for the large momentum anisotropies at early times non-perturbatively and thus allows us to start the hydrodynamic stage earlier. The viscous anisotropic formalism improves upon leading-order anisotropic hydrodynamics (see [5] for a review) by also including previously neglected residual components of the shear stress tensor through Isreal-Stewart-like perturbative transport equations. The combination of these advances leads to faster and quantitatively more reliable dynamical simulations of heavy ion collisions with fewer parameters.

2 GPU-accelerated (3+1)-dimensional second-order viscous hydrodynamics

Our implementation of second-order viscous relativistic fluid dynamics on graphics processing units (GPU-VH) is described in [3]. This code has since been adapted to evolve the second-order anisotropic hydrodynamic equations of motion, but not yet on GPUs. The performance of GPU-VH, at different spatial resolutions, is measured via the time it takes to complete 100 full time steps. Table 1 compares performance of GPU-VH on the GeForce GTX 980 Ti graphics card relative to the CPU-VH code run on the host machine with a 2.6 GHz Intel Xeon CPU E5-2697 v3 using a single core. We observe speed-up factors of ${\cal O}(100)$.

3 (3+1)-dimensional second-order anisotropic hydrodynamics

Assuming negligible net baryon density, relativistic fluid dynamics is described by the conservation laws for energy and momentum, $\partial_{\mu}T^{\mu\nu}(x){\,=\,}0$, complemented by relaxation-type evolution equations for the dissipative flows. For anisotropic systems $T^{\mu\nu}$ can be decomposed with respect to the fluid four-velocity $u^{\mu}$ and the space-like four-vector $z^{\mu}$ (defining the direction of the largest anisotropy, which for heavy-ion collisions is the beam direction [5]), parametrized by $u^{\mu}\equiv(u^{0},u^{1},u^{2},u^{3})=(u_{\tau},\vec{u}_{\perp},u_{\eta})$ and $z^{\mu}=\gamma_{z}(\tau u^{3},0,0,u_{\tau}/\tau)$, where $\gamma_{z}^{-2}\equiv 1+u_{\perp}^{2}$. Identifying the energy density ${\cal E}$ with its equilibrium form via Landau matching and demanding that the longitudinal pressure in the direction of the anisotropy is equal to its “anisotropic equilibrium” value [6], we can decompose the energy-momentum tensor as (indicating “anisotropic equilibrium” quantities with an over-hat and $\tilde{\cal O}={\cal O}-\hat{\cal O}$) [7]:

[TABLE]

Here, $\hat{\cal P}_{L}$ is the total longitudinal pressure, the transverse pressure ${\cal P}_{\perp}\equiv\hat{\cal P}_{\perp}+3\tilde{\Pi}/2$ is the sum of the “anisotropic equilibrium” pressure $\hat{\cal P}_{\perp}$ and the residual bulk viscous pressure $\tilde{\Pi}$, $\tilde{W}^{\mu}_{\perp z}\equiv-\Delta^{\mu}_{\perp,\alpha}T^{\alpha\beta}z_{\beta}$ is the energy-momentum diffusion current in the $z$ direction, and the transverse shear stress tensor is $\tilde{\pi}_{\perp}^{\mu\nu}\equiv T^{\{\mu\nu\}}$. The transverse projection tensor $\Delta_{\perp}^{\mu\nu}\equiv g^{\mu\nu}-u^{\mu}u^{\nu}+z^{\mu}z^{\nu}$ is used to project four-vectors and tensors into the space orthogonal to $u^{\mu}$ and $z^{\mu}$. By construction, the dissipative terms satisfy the constraints $u_{\mu}\tilde{W}^{\mu}_{\perp z}=z_{\mu}\tilde{W}^{\mu}_{\perp z}=u_{\mu}\tilde{\pi}_{\perp}^{\mu\nu}=z_{\mu}\tilde{\pi}_{\perp}^{\mu\nu}=g_{\mu\nu}\pi^{\mu\nu}\equiv 0$.

To close the conservation laws, additional evolution equations for $\hat{\cal P}_{L}$ and $\tilde{\Pi}$ (entering in Eq. (1) through ${\cal P}_{\perp}$), and for the residual dissipative currents $\tilde{W}^{\mu}_{\perp z}$ and $\tilde{\pi}^{\mu\nu}_{\perp}$ must be provided [7, 9]. In addition, we must determine $\hat{\cal P}_{\perp}$. In kinetic theory

[TABLE]

where the last term is the trace of the energy-momentum tensor. The $\bar{\cal R}$-function above depends only on the leading-order anisotropy of the system (measured by the macroscopic quantities ${\cal E}/\hat{\cal P}_{L}$); we take the connection between $\bar{\cal R}$ and ${\cal E}/\hat{\cal P}_{L}$ from kinetic theory. The relaxation equations can be obtained from the Boltzmann equation as [7]

[TABLE]

We refer the interested reader to Ref. [9] for a complete listing and detailed discussion of all the terms on the right hand sides in Eqs. (3)-(6). The “transport coefficients” $\bar{\zeta}_{\alpha}$, $\bar{\kappa}_{\alpha}$, $\bar{\eta}_{\alpha}$ (where $\alpha$ is shorthand for all the sub-and superscripts attached to those letters in Eqs. (3)-(6)) from kinetic theory are formally listed in Ref. [7]. Following the procedure in [8] for standard viscous hydrodynamics, we evaluate them from kinetic theory for a system of light ($m/T{\,\ll\,}1$) particles in an “anisotropic equilibrium” state, keeping only the leading terms in $m/T$. After expressing $m/T$ in terms of ${\cal E}$ (or ${\cal P}_{0}$) and $c_{s}^{2}$ we take the latter from lattice QCD [9].

To numerically solve the conservation laws $\partial_{\mu}T^{\mu\nu}{\,=\,}0$ and relaxation equations (3)-(6), we explicitly propagate $T^{\tau\mu}$, $\hat{\cal P}_{L}$, $\tilde{\Pi}$, and all components of $\tilde{W}^{\mu}_{\perp z}$ (4) and $\tilde{\pi}_{\perp}^{\mu\nu}$ (10) using a two-step Runge-Kutta scheme for the time integration and the Kurganov-Tadmor (KT) algorithm [10] for the spatial derivatives. The other necessary ingredient in the algorithm is to reconstruct the inferred variables ${\cal E}$ and $u^{\mu}$ from the numerically evolved quantities. This is more involved in anisotropic than in standard viscous hydrodynamics [9]. We first define the known quantities $M^{\mu}\equiv T^{\tau\mu}-\tilde{\pi}_{\perp}^{\tau\mu}$. Then taking the combination $(u^{3})^{2}M^{0}-u^{0}u^{3}M^{3}$ and using $u{\cdot}u=1$ and the orthogonality condition $z_{\mu}\tilde{W}^{\mu}_{\perp z}=0$, it can be shown that $z^{\mu}$ can be entirely written in terms of known quantities. Therefore, we can further define the known quantities $\bar{M}^{\mu}\equiv M^{\mu}-2\tilde{W}^{(\tau}_{\perp z}z^{\mu)}$, leading to a scalar equation for the magnitude of the transverse flow velocity:

[TABLE]

It can be solved iteratively with the help of ${\cal E}=\bar{M}^{0}-\tau\tilde{F}\bar{M}^{3}-u_{\perp}(1+u_{\perp}^{2})^{-1/2}x\bar{M}_{\perp}$, which enters in Eq. (7) through ${\cal P}_{\perp}$ and the equation of state ${\cal P}_{0}(u_{\perp})\equiv{\cal P}_{0}({\cal E}(u_{\perp}))$. $\beta_{\perp}$, $\beta_{L}$, and $\tilde{F}$ are known quantities [9]. Knowing $u_{\perp}$, all of the components of the fluid velocity can be reconstructed (with $x\equiv(1-\tilde{F}^{2})^{1/2}$) from

[TABLE]

4 Comparison

As an example, we show in Fig. 1 a comparison of the expansion velocity and macroscopic pressure anisotropy along the $x$-axis between standard and anisotropic viscous hydrodynamics for a specific shear viscosity $\eta/s=0.2$. (Many more comparison plots, including the dynamics along the beam direction, can be found in [9].) The evolution is started at $\tau_{0}=0.25\,\mathrm{fm}/c$ with a spatially fluctuating (MC-Glauber) initial energy density profile normalized such that the average over many such profiles would correspond to an initial central temperature $T_{0}=0.6$ GeV. The initial local momentum-space anisotropy is taken to be independent of position but relatively large, $\xi_{0}=10$. Fig. 1 compares ideal hydrodynamics (solid black line) with viscous hydrodynamics (green dotted line, labeled VH), leading-order anisotropic hydrodynamics (blue dashed line, labeled VH), and second-order viscous anisotropic hydrodynamics (red dashed line, labeled VAH), both at a very early proper time value of $\tau{\,=\,}0.5$ fm/$c$ (left panels) and at $\tau{\,=\,}3$ fm/$c$ (right panels). At $\tau{\,=\,}0.5$ fm/$c$ the system is dominated by longitudinal expansion, and the close agreement between AH and VAH reflects the smallness of the residual shear stress components $\tilde{W}^{\mu}_{\perp z}$ and $\tilde{\pi}_{\perp}^{\mu\nu}$. Comparing VH with VAH, we see about $15\%$ more pressure anisotropy ${\cal P}_{L}/{\cal P}_{\perp}$ in the center of the fireball and about $30\%$ larger anisotropies outside the freeze-out surface near $|x|{\,=\,}5$ fm. In fact, VH gives negative ${\cal P}_{L}$ outside the freeze-out surface which is not the case for VAH. Later, at $\tau{\,=\,}3$ fm/c, the transverse expansion has had enough time to grow to similar order as the longitudinal expansion rate, at least near the transverse edge of the fireball. At this time, the fluid velocities from VH and VAH agree well with each other, while AH is more similar to ideal hydrodynamics, indicating that the residual transverse shear viscous effects that have been ignored in AH are important for the transverse flow velocity profile. For the pressure anisotropy, the residual dissipative corrections in VAH, responsible for the difference between the dashed red (VAH) and blue (AH) lines, are much smaller than the dissipative correction in VH (responsible for the difference between VH and ideal hydrodynamics). This demonstrates the strong advantage of anisotropic over standard viscous hydrodynamic expansion schemes.